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ABSTRACT 

Digital control offers exciting possibilities for ace 
curacy and flexibility that cannot be achieved by analog 
means, The USQ-20 General Purpose Digital Computer is on 
board several ships today for the Navy Tactical Data System, 
and an effort is made to show how such a computer might alsa 
be employed directly in closed-loop shipboard Lire contrel 
systems, Some of the theory of sampled data svstems and 
digital compensation is discussed, A digital comeuter oro- 
gram which simulates the action of a sampled data system 
with a digital controller is employed to experimentally de~ 
termine the behaviour of various physical plants subject to 
Gigital contrel, An evaluation is made of the noise rejesc= 
tion capabilities, steady state error, and trensient control 
of a Gigital process, A full description and instructions 
for its use is given for the simulator program, which has 
proved to be very nelpful in making theoreticai studies of 
control system behavicur, 

The writer wishes to express his appreciaticn for the 
assistance and encouragement given him boy Professor Mitchell 
L, Cotton of the U, S, Naval Postgraduate School and JTshn 8, 
Slaughter of the U. S, Navy Electronics Laboratory in this 


investigation, 
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i °£introducwion. 

A control system which contains a digital computer 
may be classed as a sampled-data system, Such systems are 
becoming increasingly more common, both in military opers- 
tions and in the commercial field, By the use of sampling 
in control systems it is possible to build simpie, sensi- 
tive, and efficient power control devices, Sampling per- 
mits the control of tremendous power by sensitive control 
elements without excessive amplification, and aiso mini-~- 
mizes the loading effect upon sensitive instraments, 
Utilizing sampled data and digital components in a control 
system allows time sharing of important parts of the system, 
This is a very significant advantage of digital control 


ee 


systems, 
Since sampled-data can readily be coded, the Gare 
Signals in digital control systems are received and trans- 
mitted in pulse-code form, and they will provige almost 
error-free chennels for transmission through noisy media, 
Consequently, the only noise in the transmission of pulse- 
code@d signals is the quantization error, Likewise, a aigi- 
tal computer in a control system allows Gata processing of 
control information so that the flexibility and versatility 
of the digital computer can be utilized te improve the per- 
formance of the control system, Furthermore, digital cen- 
trol techniques inake feasible the system compensation by 
non-linear programming and adaptive, or self-optimizing, 


control, This is another significant advantage achievable 





weit 


Witeh digitalwcontro 
The principal features of a typical system containing 
a Gdigital computer are shown in block diagram form in Fig-~ 


ure 1, 


DG EN pee ike 
= i: 


mm T2 13 
et- 





Figure 1, Typical block diagram of a hybrid control system, 


A characteristic feature of this system is the fact that 
data will appear at one or more places within it as a se= 
quence of numbers or fixed signal levels that do not change 
appreciably except at fixed instants in time, The process 
that converts continuous data into such a sequence is called 
the sampling process, and the mechanism that implements the 
process is called an analog-to-digital converter (ADC), The 
process that converts the output number sequence of the digi- 
tal computer to a more or less continuous analog quantity 
is called the desampling process, and the mechanism accomp= 
lishing this process is called the digital-to-analog con- 
verter (pac) .727 

The quantities Tyo T5, and qT. are the periods of the 
samplers shown, In almost every case Ty) is the same as To, 
but the relation between T, and T. may be chosen by the de~ 
Signer, The simplest relation, of course, is when T. equals 


iy The system is then called single-rate since all discrete 
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number sequences appearing within the system have the same 
sampling period, A multi-rate system occurs when T, and 

T. are not the same, Analysis is simplified when Tt. and 

T are related by an integer and we have "synchronous" 
sampling, It is also simpler to analyze a multi-rate sys= 
tem wherein T3 is less than T, (ae . » T, T,/n a 
which is fortunate, since most of the advantages of a multi- 
rate system are derived when the output of the digital con= 
troller occurs at a faster rate than the input, 

The experimental results described in this thesis are 
obtained from a digital computer program that simulates a 
hybrid control system from beginning to end, See Appendix I 
for the description of this program, In an actual control 
system the only operations that would perforce be external 
to the digital computer would be the plant itself and the 
ADC, In the investigation conducted by John Slaughter at 
the U, S, Navy Slectronics Laboratory, San Diego, California, 
a Packard-Bell Model M2 analog-to-digital converter is em-= 
ployed, and a 100 micro~second conversion and settling time 
is allowed, For a discussion of a proposed analog=-to-digital 
system for use with the CDC 1604 digital computer see Ref-= 
erence 6, 

The generation of the error signal by comparison be- 
tween input and output signals can be accomplished either 
externally to the digital computer or by the digital com- 
puter itself, The latter case would require a two channel 


ADC and a multiplexing system for sampling respectively the 


UW) 





input and the output signals, 

Data reconstruction by the DAC can also be performed 
by an external device and can be of various degrees of com=- 
plexity, However, the fact that the plant in control systems 
is usually low-pass in frequency response makes the use of 
overly complex data-reconstruction systems unnecessary, 
Generally speaking, the zero-order data hold is found to be 
adequate for most systems found in practice, This being so, 
the data reconstruction process can easily be implemented 
by the digital computer since an output register of the com- 


puter can be utilized as a zero-order hold, 





2. The effect of sampling, 

Sampling is the act of examining a continuous function 
at specified increments of an independent variable, usually 
time, Sampling processes can and do assume a wide variety 
of forms and can be separated roughly into three categories: 

1. Linear, fixed pattern sampling 

2. Non-linear, signal dependent sampling 

3, Random sampling 

Included in the first category is the historically 
conventional concept of “instantaneous” sampling of fixed 
period T, By instantaneous it 1S meant that the duration, 
or pulse width, of the sample of the signal is very short 
compared to the shortest time constant of the system, Une 
der this assumption (which will be employed in this thesis) 
the sampling process can be regarded as generating a se= 
quence of impulses, This viewpoint leads to the develop- 
ment of the Z-transform calculus and also to the use of 
difference equations as a mathematical tool, Extensive 
tables of Z-transforms may be found in the literature, and 
it is relatively easy to take the Z-transform of simple, 
low-order transfer functions, However, for a fifth order 
plant or higher, finding the Z-transform is quite laborious, 
Appendix 2 describes a digital computer program written for 
transfer functions of Type O and Type 1 servo-mechanisms 
up to the tenth order, 

The Sampling Theorem by Shannon states that a band-~ 


limited signal of highest frequency Bo can be completely 





specified by sampling at a rate of 2F oe Recovery of the 
signal at this minimum theoretical sampling rate requires 
an ideal filter, In practice, signals are not ordinarily 
band-limited and filters do not have ideal attenuation 
characteristics, Even at high sampling rates, sampled~ 
data systems usually exhibit distortion in the form of 
"ripple", One means of combating this ripple is to employ 
a multi-rate digital cortroller whose output occurs at a 
faster rate than the basic error-sampling rate of the over- 
all system, By applying control signals to the plant more 
frequently much of the ripple can be reduced and the plant 
can be made to respond more guickly to an input signal, 
Design judgment must be used, however, since the output 
signals from a multi-rate controller are frequently of 
greater magnitude than those from a single-rate device, 

and the possibility of plant saturation as well as the 
added complexity of the controller must be balanced against 
the advantages to be gained from a multi-rate system, 

Of crucial importance in a sampled-data system is the 
effect of sampling on the stability of the system, In an 
otherwise stable system, sampling reduces the degree of 
stability and can even lead to instability,’-/ There are 
various criteria for continuous systems that can be applied 
to sampled-data systems, Among these are the modified 
Routh-Hurwitz criterion, the Schur-Cohn criterion, and the 
Nyquist criterion, As an example of the determination of 
stability using the modified Routh-—Hurwitz criterion we 


will use the third order Type 1 servo-mechanism shown in 
6 





Figure 2, The plant is closed-loop stable with the gain 
indicated, We desire to convert the continuous system to 

a Sampled-data system, so we first insert sampling and data 
reconstruction as shown, The matter of digital compensation 
will be held in abeyance while we investigate the effect of 


sampling only, 
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Figure 2, A third order sampled-data system, 


Suppose that a sampling rate of 100 samples per sec-~ 
ond, or T = ,01 seconds, has been chosen a priori, Then 
the modified Routh-Hurwitz procedure gives a "yes or no" 
result in regard to stability. For this value of T the 


forward path pulse transfer function (or Z-transform) is: 


P(2) 010600 2” + NOSSCGVA 2 AnewmooIws 54 





G (#) = 


a 3 
Q(z) Zee D7 SSC ea). OGD) = ON C7 
For the unity feedback configuration shown the closed loop 


transfer function has a characteristic equation given by: 


PGzm+sOC7) s=eA(z). = A383 2 tee a Al. eee 
where A3 = 1,000000 

A2 = -2,264786 

Al = 1,695799 

AO = =O 2380157 


The closed loop sampled-data system is stable if the roots 


of A(z) lie within the unit circle in the Z-plane, By 
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employing the Bilinear Transformation, Z=1 +w, the in- 
Lle=w 
terior of the Z-plane is mapped into the left half of the 


w-plane, Once the modified characteristics equation is 
obtained, namely: A(w) = B3 w? + B2 we + Bl w + BO then 
the Routh array may be formed and the criterion for a third 


order system may be applied, This criterion is: 


Bl°B2 = BO°B3 is a positive quantity if all the roots 
of A(w) are in the left half plane 


In this example, B3 = 5,340742 
B2 = 2,428516 
BL = 0,179886 
BO = 0,050856 


so B1°B2 —- BO°B3 = + 0,165247 and the system is stable for 
T equals .O1 seconds, 

If it is desired to decrease the sampling rate as much 
as possible in order to time~share the computer or to allow 
for increased computation time between samples, then it be= 
comes necessary to know the limit of stability as the samp- 
ling period T is increased, Unfortunately, the character- 
istic equation is transcendental in T, and we cannot solve 
directly for T at the stability limit, Iterative procedures 
most readily delegated to a digital computer must be em- 
ployed, Appendix 3 shows the development of the character- 
istic equation and the Fortran program for the CDC 16004 
digital computer which calculates the Routh-Hurwitz criter- 
ion for this third order system for one hundred values of T 
ranging from ,001 seconds to 0,1 seconds, From the results 
shown the stability limit on T is evidently near ,025 sec- 
onds, since this is the last value for which the R-H criter-- 


ion is a positive quantity, 
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Figure 3 shows the response of this un-compensated 
third order sampled data system to a unit step input for 
three values of T, For T of ,O0O1 seconds the response is 
indistinguishable to the eye from the results obtained in 
an analog computer simulation of the problem, This would 
be a good value to choose for the time increment ina 
solution of the plant differential equation by numerical 
methods (see Appendix 1). It can be seen that the sta- 
bility of the system is degraded somewhat when T is in= 
creased to ,0Ol seconds, but the system is still stable, 
This is the sampling period to be employed later when we 
introduce digital compensation to this third order system, 
A sampling period of .05 seconds results in severe insta-= 
bility, as we could expect from the Routh-Hurwitz deter- 
mination, Digital compensation is a must in this case, 
since stability is a prime requisite, Whether or not we 
compensate a stable system depends on how satisfactory the 


response is without compensation, 
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3, The digital process, 

As previously mentioned, a sampled-data system is 
characterized by the fact that in certain areas of the sys- 
tem signals appear as discrete quantities or fixed numbers 
in a sequence at regular intervals of time, The most cb- 
vious point at which this occurs is at the digital computer 
or controller, The input to the computer is a sequence of 
numbers and the computer output is also a sequence of num- 
bers, We can thus regard the computer as a "black box” as 


shown in Figure 4, 





E2 | eta | 
P= ee rea 


Et Lets | 
P= I; 2, ae N 






Die TA 
COMPUTER 






Figure 4, Generalized concept of a digitai computer. 


The input sequence, El [pT]. consists of N discrete 
values separated in time by Ty seconds, The output se- 
quence may contain a different number, M, of fixed quan- 
tities and the period, Ts, of the output sequence need not 
be the same as the input period, 

For simplicity let us first take T; equal to To. 
Then at a certain time, pT, the general form of iinear 
computation that the digital computer can perform may be 


expressed as: 


e2leT] + bt-e2[(e-UT] +++ Hb ed [Comet T | 


= ad-et[ptT] +at-et[(p-oT] +--+ a, pellCrnenr] 


me 





N-d et 
This can be written: e2 [ eT | = Var oa —_ 6 Ck p-; 
Wa Jif 
—K 
or, by 2Z-transform calculus: Lk (2) _ WD Zayr Zz 


4l(2#2) 4 Z bz 


which can be considered the pulse transfer function of a 
digital computer, If all the b, are zero then we have what 
is called a finite memory process since the output at pre- 
sent time, pT, is the weighted sum of the present input 
plus earlier inputs up to some finite limit, N, If not all 
of the b coefficients are zero then a weighted contribution 
of earlier outputs is also involved in forming the output 
at time pT, Since each earlier output was a combination of 
both input and output sequences, ad infinitum, we now have 
what is called an infinite memory process, Another way of 


showing this is to perform long division on the transfer 


function to obtain new weighting coefficients, Che for the 
input sequence alone, Only in the case where a2" con-= 
tains 1 ~ b,.Z J as a factor does long division fail to 


result in an infinite series, showing the need for an in- 
finite memory of previous inputs, When such a transfer 
function describes the over-all closed loop process, we 
have what is known as a "non-finite settling time" process 
and the systematic error (or difference between input and 
output) can never theoretically reach zero, 

In the system configuration employed in this thesis 


we will use the concept of the digital computer transfer 
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function in two ways, First is the over-all pulse transfer 


function of the closed-loop process, defined as; 
K(2) = aE ®, OUT FP Omi SH OUe NCE 
R(2) ~ = INPUT SEQUENCE 


Our attention will be devoted mostly to finite memory K(z) 








consisting of numerator polynomial only, Such a process 
has a finite settling time and can "forget" the effect of 
initial errors or acquisition transients in the input sig- 
nal. 

Second is the pulse transfer function of the digital 
controller in the forward path, Storage will be allocated 
for a finite number of values in both the input error se- 
quence and the output error sequence, and the controller 
transfer function will be a ratio of polynomials in z= 


defined as: D(z) = Eout(z) . processed error sequence 


Bin(z) observed error sequence 
In case the input sequence to the controller has a 
period, T 


adifferent from the period, T of the output 


iY a 
sequence we can still write a pulse transfer function that 
is amenable to interpretation using difference equations 
providing Ty and T, are related by an integer, The pulse 
transfer function is then expressed in the "Z-domain" 
created by the faster sampler, which in this thesis will 


be at the output from the controller, In the next section 


is provided the mathematical basis for these assertions, 
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4, Response criteria, 

In view of the flexibility possible with active digi- 
tal compensation there is a very large number of possibile 
over-all response functions which can be implemented, As 
a starting point, the simplest or "minimal" prototype re= 
sponse functions are convenient, Minimal prototype systems 
are approached with the point of view that they must be 
able to respond satisfactorily to some convenient test in-= 
put such as a step, ramp, or constant acceleration, or all 
three, The requirements which are set for minimal prototype 
response functions ares: 

a. The over-all response and the response of all eie- 

ments of the system must be physically realizable, 

b, The steady-state response to the test input must 

have zero systematic error, 

c, The transient response should be as fast as pos-= 

Sible and the settling time should be equal to a 
Finite number of sampling sntervals,/*/ 
The following table shows the test inputs and the required 


minimal prototype response functions, K(z): 


ere I __ 
Step aoe 

Ramp oe7t _ 37? 
Acceleration 327+ . 3277 ae 


These minimal responses are possible only for those plants 
whose pulse transfer functions, G(z), contain no zeros or 


poles outside the unit circle in the Z-plane, 
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The essence of digital compensation is that the con- 
stant coefficients of the pulse transfer function, D(z), 
for the controller are adjusted so that the desired over= 
all pulse transfer function of the process, K(z), is reai- 
ized, Depending on the structure of the sampled-data feed= 
back system, the controller transfer function is derived 
from system functions in various ways, The configuration 
employed in this thesis is the typical one shown earlier in 
Figure 1, viz., in-line digital compensation of a unity 
feed-back, error-sSsampled control system, The ADC can be 
thought of as a single sampling switch operating with per-= 
iod T, Data reconstruction will be accomplished by a zero= 
order hold, and the pulse transfer function of this DAC in 
combination with the continuous plant will be denoted by 
G(z), Since the actual input and output are continuous 
signals, fictitious samplers are shown producing R(z) and 
C(z) so that we may speak of the over-all pulse transfer 
function, K(z),. For convenience, Figure 1 is re-drawn 


embodying these ideas as Figure 5 below, 








DIGITAL PLANT CINCLUDES 
DATA D 


a ge =) iia 


Figure 5, Sampled-data system with in-line compensation, 


Taking the single-rate case first (where n=1, thus T/n=T), 
the over-all pulse transfer function is seen to be 


L® 





Clay = SD 
R(z 1 


K(z) = 
which can be solved to give 
Dele = atsy ° K(z) 
G(z l= K(z) 

Derending on the poles and zeros of the plant which 
it is desired to compensate digitally, complete freedom of 
choice of K(z) may not exist, According to Franklin and 
Ragazzini (4) the restrictions on K(z) ares: 

It is necessary that the specified over-all pulse 

transfer function, K(z), contain as its zeros all 

those zeros of the plant pulse transfer function, 

G(z) which lie on or outside the unit circle in 

the z-—plane, and that 1-K(z) contain as its zeros 

all those poles of the plant transfer function which 

lie on or outside the unit circie in the Z-plane, 

Systems which are designed for minimum and finite 
settling time often do not give good performance when sub- 
jected to an input other than that for which they are de- 
Signed, In this sense, minimal systems may be regarded 
as highly "tuned", In addition, the severe shocks which 
result in the plant cause substantial ripple in the con- 
tinuous output even though the output is error-free at the 
Sampling instants, These effects were noted in the lit-— 
Beare’ |/ and led to the introduction of a term in the 
over-all response function Known as the staleness factor, 
This factor leads to a softening of the response, with the 
result that the system can be expected to respond slightly 
Slower but adequately to a number of test inputs, The 
choice of staleness factor can be arrived at by optimizing 


Procedures” or by observing the response to some most 
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likely form of test input, The staleness factor is intro- 
duced as the constant "c" in the following expression for 
the over-all prototype transfer function: 
ee a-Klz)mininat 
(l-cz ~) 

where in practice, N = 1, ¢c is a positive constant less 
than unity, and "a" is determined by applying the Final 
Value Theorem, 

In order to compensate the system so that the output 
is ripple-free after a finite settling time the following 
/4/ 


design rules must be applied: 


a. All the rules for minimal protctype response 
systems apply. 


b, The open loop transfer function must be capable 
of generating an output that is the same as the 
input, 

c, The over-all pulse transfer function, K(z), must 
contain as its zeros all the zeros of the plant 
pulse transfer function, G(z), and not just the 
zeros of G(z) which lie outside the unit circle 
in the z-plane, 

Let us turn now to the case in which the output of 
the digital controller occurs more frequentiy than the in- 
put, The development that follows is due to Kranc (Refer-= 
ence 9), 

If we apply a discrete sequence, E(2), to a continuous 
plant, H(s), then it follows that the continuous output is: 

C(s) = H(s)-8(z) (7 
For n equal unity we may take the pulse transform of both 


sides of equation (1) for a sample period of T, Since E(z) 


is already a function of Z it is invariant under the trans- 
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formation, so the result is: C(z) = H(z)°E(z) 

When the output is sampled "n" times more frequently 
than the input (n = 2,3,...) let us adopt the notation that 
2, denotes a pulse transform taken for a sampling period of 


T/n. Then we may defines 


cz,) & Ze fete} =Z, {aed ste} 
H(Z)) Zr Hsyt 


Once again E(z) is invariant (also see Reference 4, page 88} 


Nid 


tb 


so c(Z) = H(Z_)°E(z) but we have a mixture of z and Z 


transforms, and the meaning is unclear, If we now define 
Vi STA, 
22s =e 


then it can be seen that Z@ = Z@ 


((b 


and —£(Z) = E(Z) so in the 

ti poses : ° ” n 

Z, domain" we can write c(Z,) - H(Z_) B(Z). (2) 
Returning now to the closed loop system of Figure 5 

we can see the following relationships: 


K(2)= SE and K(2n) = See 


wes 

mG) = Riz) — C@) (3) 

C (2) = E1(2) % §D(2n)S An) t (4) 

C (Zn) = E1(2)- D(2n)- (Zn) (5) 
R(z) 


Using (3) and (4): E1(z2) = (6) 


4+ KF [D(en)- Cen)5 


Substituting 
(6) into (5): _R(@)-D@n)+ En) (7) 


C (@n) = 
| 4+ 'F {D@n)- G2n)5 


18 


= = Gee se & Se 7 
SS p> ss ————EE SS Sl 


y= 
- he | 





D(Zn): C(Zn) 
ji ae Z { Dan) 6 (@n)$ 


Since K(@n)= ea there results: en) = 


which can be re-written as: 
K(2n) %, { D(2n): G(2x)} — Ol@n): © (2x) + K(2n)= O 
Equation (8) is a functional equation of the form 


Z {F(2n)} - F(2n) = — KC@x) 


that the solution is 


(8) 


It can be ar.” 


ia K (Zn) a ee 
Flaw) = SER = Dl2n) C(x) 


so finally, the transfer function of the multi-rate con- 





Broller iss - | ; K (Zn) 
V2) Fay 1 — Keo) Sal 


The ae transform of the plant is simply obtained by 
substituting 42, for Z and T/n tor FT in Gla). Hae Z.- trans- 
form of the overall process is not so simply obtained, We 
must start by imposing steady-state criteria and the re= 
straint of physical realizability. The following develop-— 


ment is due to Franklin and Ragazzini (Reference 4), 


For inputs whose Laplace transform is R(s) = = we 
S 
ite 
can wri R(2,) = P( . Aer, amd (ea) P(e) _ Cn) 
(\- 2;') (maak Kw) 
Thus the error sequence, R(Z_) - c(Z_), can be expressed as 
P(Zn) P (2) KC en) 
Ei(zn)= -—  - m ai (10) 
n) = (\-z eens (- 2n")* 
For a stable process the Final Value Theorem states that 
- lim -Z,')E1 (2 
EL(o) = sim (1-2) )E 4 (2p) (11) 


In order for equation (10) to meet the requirements of 


equation (11) we must have 
-n) IS \ 2 “N+! 
K(2@n) = ‘eye i ert tr ws ee 
\- Zy 


a2 





where f(z) is an arbitrary function of Z whose coeffic- 
ients are to be determined, 
Substituting (12) into (10) results in: 


Tee P(2n) - P(2,) + (zn) hs 


sais Nate 
(\ - B47’) 
and for zero steady state error the numerator of £1(z) 


must contain as a factor the denominator of the input, 


R(Z), Pee .’, 
Pen) - 0(20) Elen) = O-2.)“Tia@y) = CG@a ae 
or equivalently: AQ (Zn) 


—— =e for Wisi) da 
d(zy')™ l2net ay 


The restraints imposed by equation (15) are applied in order 
to determine the coefficients of the arbitrary function, 
£(Z_). Once £(2,) is known, equation (12) gives K(Z_). 

The only unknown now remaining in the fundamental 
design equation for the digital controller, equation (9), 
is K(Z,). This is obtained from K({Z_) by merely taking the 
aaa powered terms of K(Z_). eee = Neee2ho we on ae 

In very special circumstances, K(Z_) may be found in 
a Simpler way by using so-called "Single-rate" techniques 
described in Reference 4, These special circumstances are: 

a. The plant must be open-loop stable, 

oye KZ) is limited to a finite settling time type 

(i,e,, numerator polynomial only). 


ee Th 


(i) 


speed-up factor, "n", in the controller out- 
put must be selected so that all transients at 
output sampling instants are over in one period 


of the input sampler, T. This is equivalent to 
20 


Ss &-@ eae 
_—_> - oo “eo = 





requiring that "n" equal the number of terms in 
the numerator of the plant's pulse transfer func- 
tiom, G(z). 

If these circumstances obtain, and if the design goal is 

to achieve ripple-free response to a step input, then for 


a plant transfer function of 


aoa). or (b, # 0) 


it can be Pe ae that the digital controller transfer 
function is given by: 
K 
2b; Zee 
D(z) = Sa ee 
Qe (1- 2x) 
é=/ 

The main limitation of this design procedure using 
single-rate techniques is that the speed-up factor, "n", 
must be the same as the number of numerator terms in G(z), 
which may be out of the question for a high order plant, 
For the fifth order plant to be discussed later a double~ 
rate controller gave improved performance, but this con= 
troller could not be designed by the above precedure, It 
is also worth noting that the number of terms in D(z) can= 
not be less than "n", 

Up to this point the discussion of response criteria 
has been predicated upon the input being a deterministic, 


noise-free signal, but it also may be desirable to judge 


ZA 





a process on the basis of its noise rejection capabilities, 
The following remarks concerning noise rejection are due to 
Monroe (Reference 1), 

In order to consider an input signal which is contam- 
inated by noise, certain simplifying assumptions are neces- 
Sary. First, assume that the system in the absence of noise 
would perform the ideal operation on the input, ie,, the 
systematic error is zero, Second, assume the signal and 
noise are uncorrelated, Last, assume the noise is "white" 
noise over at least one sampling period, Then it can be 


shown that the mean squared error of the process is: 
N 
7 a Z 
Ks | 


where Gy is the mean squared value of the input noise 

Ak is a coefficient of the digital process, K(z). 

Ke )e= Zaz 

The ee eee reduction factor" is defined as the 

ratio of the mean squared error of the output function to 
the mean squared value of the input noise and is seen to be 
the sum of squares of the coefficients of the digital pro- 
cess, Clearly, the smaller this factor, the better the 
noise rejection capabilities of the digital program will be, 
and in general it is desired that this factor shall be less 
than unity, In passing, note that the minimal prototype 
responses to ramp and acceleration inputs have variance 
reduction factors of five and nineteen, respectively, 


The digital process can be constrained to have no 
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systematic error when the input is taken as a polynomial 
of a given order or lower, For example, if linear con= 
straints are imposed, then no systematic error exists for 
ramp or step inputs, If, subject to linear constraints, 
the variance reduction factor is minimized, the following 
expression is obtained for the coefficients of a digital 
process employing N samples of past history and precicting 


forward &€ sampling intervals, 


Ay - 2EANAzIDE EX _ E(N-1 + ax) | 
a NON +71) NON+1)(N-1) 


and the variance reduction factor is: 


S 2 2CANA1NN-1) #1IAK(N-1) # (2 x* 
ee . | NON*~7) 


The requirements that K(z) be able to predict for- 
ward by&T arises due to the delay inherent in the computing 
process, D(z). The discrete error value at sampling time nT 
equals the input signal at time nT minus the output result- 
ing from a control pulse applied to the plant one sample 
period earlier, ie., at time (n-1)T, In order to have zero 
error, the output must equal the input, so K(zZ) must precict 
forward one sample period, thus & = 1, 

The second formula above can be solved for WN, given 
any desired value of variance reduction factor, In the case 
of unity variance reduction the result shows that six samples 


are required, and the process pulse transfer function is: 
af 2 is / —é 


| —2 “3, t= mi 
K(@)s Zz > es +z + 7 = eH - z= 


a3 






eo > «¢ Baa > <=-2s «@ - - 
a alle 
=—= —-  —- Gap eae @ aay 
a Gin —_Z_Z—TT—TT  _ — a 
). SS ——S-— - 
= a. x 
@ @s @Geas=-aa- = 


__ = —_ > 





* 


~ —_— - . ee ass) 
> = «© 7 


J? Sh 


- a e. &@-%f -'* 


It was noted earlier that the presence of a plant zero 
on or outside the unit circle in the z-plane prevents the 
realization of minimal prototype response, Similarly, such 
a zero will prevent realization of the unity variance re- 
sponse, Indeed, when a zero must be accounted for, it be- 
comes quite laborious to determine any set of coefficients, 
Bs that will acrieve noise reduction, According to Monroe 
(Reference 1) it is possible to design K(z) with noise re- 
ducing properties for a plant with a single zero outside 
the unit circle by the following procedure: 


let b = the value of the zero to be accounted for, 


od 


unity for prediction forward on one sample 
period, 
The pulse transfer function of the process using N data 


points is 


N _ b, ee a 
We). 7 ae = (/- pee A? 
= / e 


K 


For linear constraints on the process (ie., zero systematic 
error for ramp or step inputs) then the matrix equation be- 
low can be solved for B A 3 a and & . These values are 


then employed in the recursion formula shown in order to 


Pan Z, and so on up Lou ame 
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It was noted earlier that the presence of a plant zero 
On or outside the unit circle in the z-plane prevents the 
realization of minimal prototype response, Similarly, such 
a zero will prevent realization of the unity variance re~- 
sponse, Indeed, when a zero must be accounted for, it be- 
comes quite laborious to determine any set of coefficients, 
Aree that will acrieve noise reduction, According to Monroe 
(Reference 1) it is possible to design K(z) with noise re- 
ducing properties for a plant with a single zero outside 


the unit circle by the following procedure: 


let b = the value of the zero to be accounted for, 


ce = unity for prediction forward on one sample 
period, 
The pulse transfer function of the process using N data 


points is Wen Z 
N kK -/ = 

Kl2z)= 2 a2 = (/- sa 
A=/ = 


For linear constraints on the process (ie., zero systematic 
error for ramp or step inputs) then the matrix equation be- 
low can be solved for 8 A 3 fa and & . These values are 


then employed in the recursion formula shown in order to 


ma rea e, and so on up to 4, 
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The recursion formula is: 
a = 
oF és i aq (bb = (iy, = 
6, = - Gilg . ecb -S4 - Fr 


Even with a noisy input signal it is still desirable 
to exercise some control over the transient behaviour of 
the system, It is possible, again according to Monroe, to 
design the over-all process according to a composite cri- 
terion wherein preference can be expressed for step tran- 
Sient response (in the mean squared error sense) over noise 
rejection by means of an arbitrary parameter, 4, ranging 
between zero and one, Subject to linear constraints and 
for no plant zero outside the unit circle, the process co- 


efficients, a are calculated from the following formulae: 


Kk? 
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Ciera, = a + 14:., aa A; (recursion relationship) 


A é. 20m Be 4 pe Ey = Ny Vo ( starting 
44-4. 44% — 4 Formulae ) 
ZB = ee | ye) = ie —/ 
oO — LO 
/ at 





by, = Go) Ne N do Lx ON. )) = ANN Ny 
=" (A, -1)* 


xf += a £ faba htf4 + ” Nigra A, f3 
wii ae 
. {2 ‘aha 7 sats) 4 +45 }* 


Figure 6 shows the Fortran coding of these equations 


Ns 


and the results for a ten point process with beta of 0,75, 
The program may be used for another design merely by chang= 
ing the two statements giving the values of the parameters 
beta and AN, 

For the design of a process to a composite criterion 
when the plant has a zero outside the unit circle see Ref- 


erence l, 
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Program for computing composite coefficients 


| Figure 6. 





5. Program requirements, 

Before discussing the realization of the digital com- 
pensator, D(z), by programming of a digital computer, iet 
it be noted that the pulse transfer function of the compen- 
sator can be realized in at least two other ways if a digital 
computer is not available or economically iiectinie,/ °7 First 
is the realization by a delay-line network, sometimes refer- 
red to as a pulsed~data processing unit, through the use of 
delay elements, potentiometers, and summing amplifiers, The 
pulse transfer function can also be realized by a pulsed- 
data network consisting of resistors and capacitors, In= 
Guctors are avoided, because in the usual frequency range of 
Gigital and sampled-data control systems the sizes and weights 
of the reguired inductors would be prohibitive and the losses 
would be excessive, Pulsed data RC networks have the advan-= 
tages of simplicity in structure and economy in implementation, 

When a digital computer is available, realization is 
relatively simple by programming, The program can generally 
be carried out by means of one of three GQifferent methods: 

a. Direct programming 

b, Iterative, or cascade, programming 

cc. Parallel programming 

Let the pulse transfer function (or Z-transform) of 


the digital controller be 





a! _kK 
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a 
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Bae first tem insehe Geneminatoreof Déz) must be non-= 

zero constant for the sake of physical realizability, In 
the direct method of programming it is necessary that this 
term be unity, since the direct method employs the differ- 
ence equation equivalent of the above pulse transfer func~ 


tion, Writing the difference equation we get 
M-! 


N= 
Eout [nt] = DAw Ew [(n-)T] — 2 B; our [r-j)T] 


i 

In the digital computer this method requires 2Mt2N-1 
cells of storage for data and constants, MtN-l multiplica-= 
tions, M+tN-2 additions and/or subtractions, and MtN-2 data 
transfers, 

The Cascade method is based on the fact that D(z) is 
factorable into its poles and zeros such that a cascade of 
first order transfer function can be formed, In the computer 
this method requires 2M+*N+3 cells of total storage, M+tN+1 
multiplications, M+tN additions and/or subtractions, and 
Mt2 data transfers, 

The parallel method of programming is based on the 
partial fraction expansion of D(z) and requires 3M+2 cells 
of storage, 2M multiplications, 2M-l1 additions and/or sub-= 
tractions, and Mt2 data Sraitae = Aue 

Before a choice of methods can be made on the basis 
of storage requirements or computation time, D(z) must be 
known, Since any advantage of one method over the others 
depends primarily on the number of terms in the numerator 


of D(z). Generally, other considerations prevail. For 
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example, in the cascade method each stored constant is either 
a pole or a zero of D(z), and thus this method lends itself 
readily to experimental manipulation of the pulse transfer 
function ina "cut and try" sense, If storage space is at 

a premium, as in a special purpose digital computer, the 
direct method enjoys the advantage that data transfers can 
be effected by a single command that processes all the data 


“7, 


in a shift register, In the simulation program described 
in Appendix 1 the direct method of programming the digital 


controller is used, 
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6, Implementation of a double-rate digital controller by 
direct programming methods, 

Direct programming takes the difference equation 
equivalent of the pulse transfer function and interprets 
a to be a delay of one computation period, Thus, in the 
Single-rate controller both input and output discrete se- 
quences are equally spaced in time, and the programming is 
straight-forward, The computer program for a double-rate 
controller will be described in terms of the various arith- 
metical and shifting operations that must take place and 
not in terms of any particular programming language, The 
method will be shown using an example, and can be easily 
extended to the case of an n-rate controller, 

Suppose that the pulse transfer function of the double- 
rate controller is: 
Dl@)= AL4 Ae A3-@," + Ad: Ba + AS. a + Aerial 

1+ 82.2; + B32," + B42, +85 2, + BezZ; 
The output from the controller will be a sequence, zout(Z,), 
of discrete values separated in time by T/2 seconds, The 
input sequence is spaced at intervals of T as determined by 
the error sampling rate. When we write the difference equa= 
tion to find a new value for Eout by cross-multiplying the 
transfer function and treating ia as a Gelay of T/2 seconds 
we get: 
Eout(n-%] = At: EN(N:% [+ AL'Em LN-)%] + AS Ew lOr-2)% | 
+ A4- Ein((n-3)% | + AS-Ewl(n-4)% | + AG- En lw-s) %, | 
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— { 62-Eour (cn-1) % | + B3- Eour [(n-2)% | + B4 Eour [(n-3)% | 


+BS: Eout[(n-4)% | + BG: Bout] (N-5)% | : 


Since the discrete sequence, Ein(z), is at intervals of T 
seconds only half of the terms involving Ein above will be 
non-Zero, Which half depends on whether N is even or odd, 
When N is even (ie., an output value is generated at the 
"same" time the error is sampled) then the terms involving 
Al, A3, and A5 are non-Zero, When N is odd (ie,, an output 
value is produced by the controller between input samples) 
then the only terms with non-zero values of Ein are those 
involving A2, A4, and A6, This amounts to alternately 
multiplying the input sequence by the even numbered coef= 
ficients of the numerator of D(Z.), then the odd numbered 
ones, ..,etc, 

In a given period of time the output sequence from a 
double-rate digital controller will consist of twice as 
many discrete values as the input sequence, For the sake 
of illustration let us number both sequences backward in 
time in the same way, thus; The most recent vaiue is num- 
ber one, the next most recent value is number two, etc, 
The total number of values in each sequence that must be 
stored in the computer memory is determined by the number 
of terms in the controller transfer function, In the ex- 
ample of this section we require Output values numbers one 
through six and input values one, two and three, Stach out- 


put value is employed only once in each numbered position 
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and passes from storage after six intervals of T/2, Each 
input value is employed twice in each numbered position; 
once as a multiplier of A; ana T/2 seconds later as a mul= 


eiplier of A Since the numerator of D(z.) contains 


isis 
six terms, each input value also passes from storage after 
six intervals of T/2, A simple-minded pictorial represen- 
tation may serve to clarify these ideas, Let the present 
time be the 16th instant of the closing of the output samp- 
ling switch of the double-rate controller, Neglecting 
computation delay (which is usually at least an order of 
magnitude less than the sampling period) the present time 
is also the eighth closing of the input sampling switch, 
The input and output sequences may be thought of as exist-= 


ing in storage as shown in Figure 7, 


input sequence, £in 


observation 


“4 


8T le CT time 
to be computed Output seguence, Eout 





present (older values) 
time 


Figure 7, Input and output sequences for D(Z,) 


In the difference equation the terms involving the input 
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sequence would be (since in this example, N = 16): 
Al-Ein(8T) or Al-Ein(1) 
A3e¢Bin(7T) or A3-Ein(2) 
A5¢Bin(6T) or A5-Ein(3) 
The value of Ein(1l) has just been measured by the input 
sampler and the controller produces Eout(1). 

Now let us move on in time by T/2 seconds, It is 
necessary to re-number the output sequence since a new 
Bout(1) will be produced, The previous Eout(1) becomes the 
current Eout(2), etc,, with the previous Eout(6) being dis-— 
carded, We are now between input samples and no rew input 
information is to be received, The previous input values 
have become T/2 seconds older, however, and the sequences 


in storage might be pictured as shown in Figure 8, 


Input sequence, Ein 





ST 7T 2 ~ 6 time 


Output sequence, Eout 
to be comouted 





present (older values) 
time 


Figure 8, Sequences between input samples 


Since the difference equation is in terms of T/2 the non-zero 
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terms involving the input sequence have changed and are 
meow (for NM ="17): 
A2°¢Ein(8T) or A2¢fin(1) 
AA@Bin(7T)oor AdeSin(2) 
A6*Ein(6T) or A6-Ein(3) 
The controller computes a new output value using all 
Se the denominator coefficients of its transfer function 
put only the even-numbered numerator coefficients, As 
before, the output sequence must be shifted since we are 
only half a period away from receiving a new observed value 
Beomuebin(l1). 
For an example of the implementation of these ideas 


using the CDC Fortran programming language see Appendix 1, 
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ome DoLications, 

The results to be described in this section were 
obtained using the digital computer program simulation of 
a control system described in Appendix 1, The restonse 
curves were obtained through use of a Graphplot computer 
routine for the CDC 160 digital computer devised by Lt, 
Robert L, Hogg of the U, S, Naval Postgraduate School, 
Monterey, California, The program in Appendix 1 also pro- 
duces numerical results on a printer, but these are not ine 
cluded in this work, 

Initial investigation was made of a simple plant 
Gescribed by a second order differential equation whose 
pulse transfer function had no zero outside the unit circle 
in the 2-plane, Next, a third order plant with a pulse 
transfer function containing a zero outside the unit cir- 
cle was investigated, Both of these plants were theoretical 
and in neither case was an actual physical counter-part at 
hand, 

At the U. S, Navy Electronics Laboratory, San Diego, 
California, a group headed by John B, Slaughter is engaged 
in studies of Digital Sampled-Data Feedback Systems for 
Shipboard Control Applications (NEL problem N4-17). The 
author was privileged to spend ten weeks with that group 
early in 1963, during which time a twin 3"/50 gun mount 
was in the process of being installed to serve as the plant 
in an experimental control system, A USQ-20 general purpose 


digital computer is to be employed as the controller, The 
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manufacturers of the gun mount supplied a fifth order 
transfer function for the gun in train, The majority of 
the simulation studies to be described in this thesis were 
cenducted using this transfer function, and it is felt that 
the results obtained are applicable not only to this par- 
ticular obsolescent plant but also to present day ship= 
board weapons, especially missile launchers, 

The simulation program has the capability of testing 
a system for input signals of either unit step, ramp, or 
acceleration, In addition, additive Gaussian noise of 
selected mean and variance may be included in the input, 
The first action of the progrem is to produce the response 
obtained in ent case of the continuous, un=-compensated sys= 
tem for purposes of comparison, In the event this response 
is unstable no graph is made of it, Instead, the test input 
is reproduced, The next action of the program is to simulate 
the same plant incorporated in a sampled-data system with a 
Single-rate digital controller and produce the response for 
this case, The orogram is further capable of simulating a 
system using a double-rate digital controller to produce a 
third response curve, although this is not done in every ex- 
periment, For convenience, the program will also compute 
the pulse transfer function of the controller when requested 
to do so, or the controller coefficients may be provided by 
the user, The Z-transforms employed in this work were com- 
puted by the program described in Appendix 2 and were veri- 


fiei by hand for the second and third order systems, 
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The second order plant has a Laplace transfer func- 
tion given by J 
SOs cS) 
The pulse-transfer function of the plant and zero-order hold 


combination for a sampling period of T = ,O1 seconds is 


= 


Bee) = Sake (/+ 0.9204552'') (OSI ee ned Coes 
(1-2°')(/- 0.778801 27') 1 = 1.772801 B7'+ 0.778 E0127! 


Mivemwe are no zeros of G(z) on or outside the unit cirele in 
the Z~plane, so minimal prototype compensation can be a- 


chieved, As noted earlier, the minimal prototype over-all 


pulse transter function for a step input is K(z) = Ja sO 
D(z)= ne Kv (z) j= pa7neeovel'+ .77 2g0n2 2 Sea 
Zz a—————— SO a So ”~—~C a ne ee ie 
Giz) |-K(2) OS Bea!  FOm Page a /~ 27! 
/~, 7788 0/27" ae _ £8.93 9394 — 14, 7500/9 27 ' 
a. qloassa-) Trait OO & O.9204 5527! 


(note that 1-K(z) contains the pole of G(z) lying on the 
unit circle, ) 


Fisure 9 shows the plant behaviour in a continuous, un-= 
compensated system and the behaviour in a sampled-data 
system using the digital controller above. Note that there 
is zero systematic error at the first sampling instant after 
the input was applied, and also at all succeeding sample 
instants, For a K(z) containing "n" terms there will be 
zero error for a design input after "n" sampling instants, 
The plant behaviour between samples is another matter, and 
in this example probably would not be considered satisfac- 


tory, unless the only desire was to increase the speed of 
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response with no regard for over-shoot and inter~sample 
ripple, 

Simulated Gaussian noise with a standard deviation 
of 0,204 was added to the input and the experiment was re- 
peated, Figure 10(a) shows the results, According to Mon- 
roe the process should have unity variance reduction factor 
Since the sum of the squares of the process coefficients is 
one, Judging from Figure 10(b), which shows the noisy test 
input, this is approximately correct, since no noise reduc~ 
tien 1S apparent, 

In an effort to improve the transient behaviour for a 
step input a staleness factor of 0,5 was selected, The 


over-all process was designed from 


al Kl@) minimal 
K(#)= 





/-0.5 27' 
and the Final Value Theorem showed that for zero steady- 
O.s2' 


Soccer a=] = c = CyS. somthat KGi= [22s 
{/-O.sz27! 


z= Osz2'! 
(1- 2'(/- 0.788012) Fog eam 


OSied | ote OGD Ct! a= Oe Sms 
J-O.s27!' 


Seay = 





/ 


= OS-0.39940052 9469697 —7,37500952 
ORE /-. ORE? POOP OTF LOS 5 sae 








The step response using this digital controller is shown 
in Figure 11, This is an example of a non-finite settling 
time process, but for practical purposes Zero error is 
achieved at about the seventh sampling instant and there- 
after, and the inter-sample ripple and over-shoot is much 
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reduced compared to the minimal design, 

This same staleness factor design was subjected to 
a noisy step input with the results shown in Figure 12, 
Since K(Z) can be expressed as an infinite series by long 


division, thus: 0.5271 + ,2527 Sas nee 


ee eS 
the variance reduction factor is approximateiy one-third, 

It is not apparent from the simulation results that the 
process indeed had this noise reducing property, 

The six point finite memory process for unity variance 
reduction (see page 23 above) was then subjected to a step 
input with and without noise, Figures 13 and 14 show the 
results, This process and the following one depend on the 
pulse transfer function of the plant not having 4 zero out= 
side the unit circle, It can be seen in Figure 13 that the 
transient behaviour of a design whose sole criterion is 
variance reduction is rather poor, For this reason the next 
process considered was a composite one in which three-quarters 
of the design effort was towards transient control and only 
one-quarter towards noise reduction, The program on page 
27 was used to compute the process coefficients, as shown 
on the same page, This gave the most complex controller so 
far with ll numerator terms and ll denominator terms, The 
results of the simulation experiment with this controller 
are shown in Figures 15 and 16, The transient behaviour 
displayed in Figure 15 is a noticeable improvement over that 
of Figure 13 where the only aim was noise reduction, 


At this point a certain disenchantment with the notion 
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of noise reduction begins to set in, The plant is a simple, 
well-behaved one, and in the four experiments so far using 
"clean" input signals both the continuous system responses 
and the sampled-data system responses can be verified ana= 
lytically with some effort, When the unit step input has 
"white" noise added to it as shown in Figure 10(b) two ten= 
tative conclusions suggest themselves, First, the contin-= 
uous system with no compensation is relatively unaffected 

by the presence of the noise, A close inspection of Figure 
10(a) reveals that the peak overshoot still occurs at about 
one-tenth of a second and its value is increased only slight= 
ly (from 29% to 37%). A steady state error equal to the mean 
squared value of the noise can be inferred from Figures 12 
and 16 in which the time scale has been doubled for a longer 
look at the plant’s response, Over=all, however, the low= 
pass character of the plant effectively prevents any drastic 
reaction to the presence of the noise, Not so in the case of 
the sampled-data system, so our second thought at this point 
is that the insertion of the digital controller in the for-= 
ward path disrupts the low-pass filtering action of the plant, 
The noise in the error signal is the same as in the input, and 
the noise power in each discrete error sample is integrated 
over the next sampling period and applied to the plant, This 
integrated noise power causes some fairly drastic gyrations 
in the plant regardless of the design criteria, A brief com- 
parison of the four noisy experiments up to this point provides 


small reason for employing a variance reduction design in 
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preference to a noise-free criterion, It appears that the 
effort could better be expended in pre-filtering or other 
measures to eliminate noise prior to the input to a sampled- 
data system, 

The next four simulation studies were conducted using 


a third order plant whose transfer function is given by 


H@® = TOROS 


SCS 255770) 


For T = ,01 seconds the z-transform of the plant~and-hold is 


7 S & 
0/06002 #¢. 033672 2 * y+ .006s e427 


G(#) = ‘a o:  -- |, ==, nr 
emer AIS SEG 2 Ff) CO6A/2 7 2 “Oreo Agee 


Gr in factorea form 


Kle)= 10lOS 2/2 9. 2073272" / + 2. RAG) ) 
(1-2-7) (/ - 0.772 2012 -"N/ — 0.4965 9S 2~') 

The presence of a zero of G(z) outside the unit circle is 
noted, This zero will prevent us from realizing a "minimal" 
prototype response, In accordance with the rules given on 


page 16 the design of a process giving prototype response 
to a step input proceeds as follows: 
Kl2)= (l#2.967A75 27) La 2 
I-K(@2) = (1- 2°) /+ b,27') 
Solving for the two unknowns, ay and by, results in 


Ge .ASAOCL Ey =. IAZ7S © and thus 


a 


/ 


P= OA9S5206282°' + 0.747938 2 


= (/ + 2. 96747527 ')0,a5 a0622~' 
Solving for the digital filter coefficients by hand, we write 
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CS, /-,778 Soa U/-.d 9605252 |) ([naGeraice ). AS LOG 2277 


D(z) = (O/06N/ +, 2093219 Le 2P6PTTSE ) (BL) / +. 74793¢2°') 


and after cancellations of the factors shown, we get 


D(a) = 22:779434 ~ 30.327923 2 "+ 7/76 S592 * 


BORD ONGS Tie +f O. /Samnoees = 
The un-compensated response and the response of the sampled- 
data system using the D(z) above are shown in Figure 17, 
The non-minimal nature of the digitally~controlled response 
is characterized by the fact that it took two sampling per- 
iods (.02 seconds) for the output to reach zero error at 
sampling instants, There is a slight amount of inter- 
sample ripple, but the peak overshoot is greatly reduced 
from what it was for the continuous system, The improved 
damping of this third order sampled-data system compared to 
the second order system considered previously is attributed 
to additional filtering action of the pole at -70, Since 
this design was for a step input there will not be Zero sys-~ 
Eemacac error for higher order inputs, such as a ramp. this 
can be observed in Figure 18 where the curve labelled Dl is 
the response of the digitally controlled system to a ramp 
input, (Note: in this figure, and all others for which the 
test inout is a ramp, a 45 degree slope from the origin is 
the input signal), The controlled response settles down to 
a steady state error of about ,017 units after two sample 
periods, Whether or not this is an intolerable error in re- 
sponse to a velocity input depends on the system, For ship-~ 
Beema £ire—control application such as tracking a target, a 
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design for zero error to a ramp is usually preferable, 
More about this later, 

An additional refinement of the step response of this 
third order system can be achieved at the cost of one extra 
period, T, in reaching steady state by designing for ripple- 


ree response as discussed on page 17, For this case we get 


Kl2) = (1 + 2.209327 BV + 2.967275 2 ')a Zz 


/- KC2) = (/- = Ys ~ O,2~° + bz any 


k 


and by equating coefficients of like powers of z™~ (k=1,2,3) 


we obtain the below equations to solve for al, bl, and b2, 
Gap = i - b, 
O.€2//3208% &, = b. 


SsOlving these, we get al = .2084318 so that the process is 


1 4 + 129463427 


imez ye 20643182  ° &4+«+ ,662105z 
These process coefficients as well as the pulse trans-~ 
fer function of the plant were provided to the simulator 
program which then calculated the digital controller and 
produced the ripple-free response shown in Figure 19, 
Figure 20 shows the same system with a velocity input, The 
steady state ramp error has increased slightly, from .017 to 
~.020 due to the additional term in the process z-transform, 
and the controller had a total of six terms, In the absence 
of input noise we have demonstrated satisfactory digital con- 


trol of this lightly damped plant when the input is a unit 


step, 
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Turning now to the twin 3"/50 gun mount, the transfer 
function of the forward path train drive is given by 


Hs) = |.606423Ex/0" 
SOs 4riXs+/7Xs+23X 5% 23) 


Most of this very large forward gain 1s suppiied by the 
amplidyne drive units and results in a quickiv responsive 
system when suitable compensation is added, In actual use 
the train drive function was elaborately compensated by 
various continuous networks in a multlLtude of feedback 
paths in order to obtain stable operation, All of this 
continuous compensation has been removed Erom the gun by 
the group at NEL, and the necessary control wiil be achieved 
by digital methods, Figures 21 and 22 were provided by 
Donald Lackowski, an engineer in that NEL group, and show 
the root locus and Routh criterion analysis For the forward 
transfer function with unity feedback, It ean be seen that 
the stability limit on the gain, K, is reoughiy three orders 
of magnitude lower than the actual value present, and thus 
the plant is highly unstable without compensation, This is 
borne out by the simulation studies to be described next, 
In every case the continuous response went rapidly cff scale 
and could not be shown in the figures to follow, As previ- 
ously mentioned, in this event the simulator program re- 
produced the test input in lieu of the unstable response, 

The computer program described in Appendix 2 was ex-= 
ceedingly helpful at this point since finding the 2-trans- 
form of this fifth order equation for just a single value 
of T is quite laborious by hand, With the computer a 
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number of transforms was quickly obtained for a spectrum 
of sampling rates ranging from one hundred samples per 
second to as slow as one sample every four seconds, For 
the faster rates the pulse transfer functions themselves 
were Of the fifth order (ie., five numerator terms and six 
denominator terms), As the rate was decreased it was ob- 
served that the higher ordered terms in both numerator and 
denominator of G(z) became progressively smaller, For T 
of one-half second or larger a gocd approximation of G(z) 
can be made usSing only a second order function, In addi-= 
tion, the zero of G(z) which lay outside the unit circle 
for smaller values of T had moved inside the circle for 
T= 0.5 seconds, The computer caiculated pulse transfer 
function of the plant and hold combination for T of one= 


half second is given below: 


Numerator of G(z) Denominator of Giz) 
1336, 331533192 1, 00000000 
897.96714053z"° ~1,00430037z7+ 
6,02702748z > 0043012527 ° 
0004922027" -,00000087z" 7 
,00000000z7> _0000000027* 


000000002"? 


The computer program also factored the numerator and gave 
the following real values for the zeros of G(z)s: 
-,Q00000001, -,00008267, =,00669675, -.660518491, 0,0 

It was decided to investigate the advantages of a multi- 


rate sampled data system, the analysis and design of which 
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was discussed beginning on page 17, The simulator pro- 
gram has the ability not only to produce results for the 
continuous system and a single-rate digitally controlled 
system, but also, upon request, to simulate a system using 
a digital controller whose output to the plant occurs 
twice as frequently as the input error sampies are sup= 
plied to the controller, This is calied a double-rate 
controller, The output rate was chasen to be two samples 
per second in order to take advantage of the second order 
approximation to G(z) mentioned above, It was desired 
that the over-all process have zero steady state error to 
a velocity (ramp) input, Following the method cf page 19 


the design for prototype ramp response proceeds as [OLiows 


The fundamental design equation for the double-rate con]= 





troller is K (22) 
Di(Zelt) nx .—X—— ae 
Elen 1S xe 
—2 
For T/2 = 0.5 seconds, G(2,) ~ /336 zZ3' + &7& Za 
/— 39 
The Laplace transform of the input is R(s)= I , whence 
Sr 
—| 
rh Ha a 2 
Ri@z)= 2°F2 5 R(zz)= #2 = R(H) 


Applying steady state and realizability constraints, we get 


Pe ™ a 
K(@a)= Qo 2a) Eye) wrote Flea)= Qrta'+ Be 


and a look ahead at the constraint equation, 


d“Q @a) 
d (z,)™ 


a C) for rsa eek ot 
a; ou 
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shows the need for only two free constants, ay and Aa59 SO 


Q(22) = Pl2,)- Pl) (2,) =] 2;'- 1-2; (ay27' + Gaz?) 


0.52, - 277 (A425' + G25") == 0,52, - ley? *- biove 2 


for Mm O: Q(z) = 0.5 —G,—a, = ae 


1= 





for n=i1 4Q (2a) = 0.5 —- 34, — 4 oye 
ee Zo=1 


solving these equations we gets 





a, = Loe 


so that K(25} = (1228 \"(\.s23'-2,") = (1+23')"(.seze Se 
Te) 


= /. § B;'+ 2.02,*— O52, — fe Zz," 





To find K (23) we take the even~powered terms with the re- 
Sunita K (22) eC) a AO) i = K (=, ) 
which can be recognized as the minimal prototype ramp 
response function shown on page 14 expressed in the Z. 
domain, 

Substituting now in the fundamental design equation, 


the pulse transfer function of the double-rate controller 


(j- 22") } F 25 » HOS, COG es Oe) 
oa = a 


(13362, ‘+ 0982,) (elias A.0O2777* + 025 sy 


Ve 
1336 - 438 27' — 298 25? 


We must normalize this result with respect to the leading 
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denominator coefficient, since direct programming requires 


eevelue Of unity in this position, The final result is 


Diz.) = 2022754 ~ .000742S50 Zo. 
1.0 — 0.327 8443125 — 0.6721ISS69 Z,* 





Before proceeding with the simulation run we must 
provide for the single-rate action with T = 1,0 seconds, 
This is simply done by providing the simulator program with 
the ramp minimal prototype function and G(z) coefficients, 
For this sampling rate an even better second order approx= 


imation is available for G(z), namely: 
+ Jae = 


een 


~ 2 
G(#)= - 

Tne responses of the two systems for a noise-free 
velocity input can be seen in Figure 23, The single-rate 
controller gave zero error at the second sampling instant 
and only moderate overshoot and inter-sample ripple, The 
double-rate system was 0,5 seconds faster but over-shoot 
Was worse due to the larger control forces involved, The 
ripple in the double-rate case was more persistent, also. 

These same systems were then subjected to step inputs 
with the results displayed in Figure 24, The highly tuned 
nature of any prototype design is exhibited here, where 
the input is of lower order than that for which the system 
was designed, An interesting result in Figure 24 which 
was not expected is the slightly decreased peak over-= 
shoot in the case of the double-rate controlier, 

The next design criterion applied was ripple-free 
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step response using both single and double-rate control- 
lers and the same error-~sampling rate of the preceding 
design, The riople-free process, K(z), was computed and 
supplied to the simulator, which calculated the single- 
rate controller, Dl. The double-rate controller was de- 
sicnec using the techniques discussed beginning at the 
bottom of page 20, Both controllers were quite simple 


and are given below; 


Dy = 200022222 2 py _ .0004475 — .00000/9aS5 #, 
[ao 2osz3a32' ° 7 /.0 


Figure 25 shows both systems responding without ripple to 

a noiseless step input, Again, the double~rate system is 
about one-half a second faster, The ramp response is next 
shown in Figure 26, Both systems have a steady state error 
Since they were designed for a step input only, in the 
case of the double=-rate controller the sustained ripcie 

is due to the fact that the design technique was a special 
One, explicitly for ripole-free step response, and cannot 

be expected to give conventional behaviour for a ditferent 
input, 

The results of the simulation studies using the seconc 
order approximation designs with inputs contaminated with 
Gaussian noise tend to reinforce the tentative concliusicns 
mentioned earlier, It was impossible to observe the response 
of the continuous system to a noisy input since the contin- 
uous system was unstable, The disturbing result is that 


the sample-data systems, both with single and doubie-rate 





Gigital controllers, which had given stable performance 
of varying degrees of goodness when the input was nolise= 
free were, with but a single exception, unstable in the 
presence of noise. The exception is shown in Figure 27 
in which the design of Dl and D2 was for ripple-free re- 
sponse to a step and the input was a noisy step, It seems 
apparent at this point that a large signal-to-noise ratio 
is a sine qua non in the input to a sampled-data system, 
In the preceding experiments it was only the pulse 
transfer function of the plant and hold combination that 
was approximated in order to simplify the design caicula~ 
tions, The differential equation of the plant itseif 
remained of the fifth order for the numerical methods 


solution of all response values, 
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The simplicity of the design techniques for sampling rates 
of one or two samples per second which result in valid 
simplifying approximations is quite appealing, but there 
lurks an intuitive susoicion that these rates are too slow 
for a system required to accurately follow a high-speed 
aerial target. This suspicion has been confirmed in pri- 
vate communications between John 3, Slaughter and Gene F, 
Franklin, the results of which are that the damped natural 
frequency of the gun mount in train is about one and one-half 
cycles per second and that sampling should occur between six 
and twelve times the damped natural frequency of any system 
for good results, For these reasons a sempling rate of ten 
samples per second has been decided on by the group at N&lL, 
and the remaining simulation studies of this work wili use 
this rate, 

For it = 0,1 seconds the pulse transfer function of 


the plant and zero-order nold is shown below: 





Numerator of Glz2) Denominator of G(z) 
21.9124153127+ 1, 00000000 
1 speenecess2- Seoueos a7 me 
63.9418040127? 2a7eses- 
3,0843658927" Sais Fera2 ° 
5 é 


, 004190482 — ,006126222 — 
~,0000015227” 

The roots of the numerator are located at the origin and at 
-,00139908, -.05255621, -.47316025, and -5,4S614151 


The presence of the zero of G(z) outside the unit 


Ad 





circle complicates the design effort slightly as it did 
in the third order plant considered earlier, The same 

methods apply and the resulting pulse transfer function 
of the process which gives prototype response to a step 


input is K(z) = .154174866272 ° 


+ ,845825134z 
Figure 28 shows this process responding to a noise- 
free step input, and Figure 29 shows the same process with 
a noise-free velocity input, The remarks made earlier con- 
cerning prototype processes apply and will not be repeated 
here, Figure 30 shows the results obtained with this de- 
Sign when the step input contains additive noise, The 
best that can be said about the plant's behaviour in this 
case is that it remained stable and in the vicinity of the 
input signal. Not so for a noisy ramp input, to which the 
response of the sampled data system was unstable, It ar- 
pears that the sensitivity of a sampled data system to 
noise increases with the order of the input signal, 
The prototype design for Zero steady state error in 
response to a velocity input begins with the equations below: 


Mea = Gee 5.48614152°*)(a,z~ 


1 - K(z) = (2 - 271)%(1 + by27) 
From these are obtained three equations to be solved for 
three unknowns, namely 


=2,=—- 5 -5,4861415a. = bio and 


i 1! 
a5 = 9.486141 5a, = 2b, -1 
Solving these, the resulting over-all process is 
1 2 3 


K(z) = .438754756z. + 2,1224904882 “ - 1,5612452442 
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The coefficients of this process as well as the coeffi- 
cients of the numerator and denominator of G(z) and the 
value of the zero outside the unit circle were supplied 
to the simulator program, which then calculated a single- 
rate digital controller, The responses of the system to 
a noise-free ramp and a noise-free step are shown in Fig- 
ures 31 and 32 respectively. With the addition of noise, 
instability again resulted when the input was a ramp, while 
for a noisy step the output was at least bounded for the 
duration of the experiment, although it began to look in- 
auspicious towards the end, See Figure 33, In an attempt 
to achieve at least a semblance of stability in response 
to a noisy ramp signal the noise deviation was reduced 
from 0,2 to .05, but to no avail, A future investigation 
will be made into the stability threshold for noisy ramp 
signals, 

It was then decided to design for prototype step re- 
sponse using a Staleness factor of 0,5, The design steps 
are skown below: 


1 + (845825134z27° as found 


Prototype K(z) = .1541748662 
earlier, This function contains the zero of G(zZ) located 
outside the unit circle at zZ = 5,48614151, The process 


pulse transfer function with a 0,5 staleness factor is 


given by p(z) = B(O./S4174 266 27'+ 0. 845895134 27) 


J/- OSza~’ 
For a unit step input, application of the Final Value Theo- 


rem shows that "a" must have the value 0,5, Thus K(z) is: 


Ts 





/ 


(077087433 27 + 0.422914567 ZB 


c=) = =, a? 
. (2), KC) 
Substituting in D(z) = oS I\- K C2) 


where P(z) is the numerator of G(z) given on page 71, and 


Q(z) is the denominator of G(z) given on page 71, 


and cancelling the common factors (1-27+)(14+5, 48614151271) 


a 


the single-rate digital controller is given by 


Numerator of D(z) Denominator of D(z) 
003517980 1,.00000000 

-,002167290271 95013414271 
00039630727 24861928277 

~,00002151027 0108837027? 
000000025274 0000103927" 


This pulse transfer function was supplied to the simu- 
lator program; test inputs of noise-free step, noise-free 
ramp, noise-free acceleration, and noisy step were applied 
with the results shown in Figures 34 through 37 inclusive, 
The main point of interest in these results is that step 
response in the presence of noise for this staleness factor 
design is perhaps the best-behaved so far, although it is 
not good by any means, In addition, noisy inputs of ramp 
and acceleration were applied, but both of these gave un-| 
stable responses as in previous designs, 

The final design effort to be described in this work 
Was directed toward achieving ripple-~free response to a 
velocity input when sampling at a rate of ten samples per 
second, For this design all of the zeros of G(z) must be 


contained in K(z), so we begin by writing 
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ty (1+, 4731662527") (1+.0526 562120) 


~2) 


ez) = (1+75,48614151z 


(1+, 001399082) (a,z7} ee 


2 
Since we are designing for zero steady state error for a 
ramp input, then it must be true that 
-1,2 -l] —2 oa 
oa 
“C1 biz + bo2 + bz 
After multiplying out and equating like coefficients of z 


game(z) = (l-z 3 


we obtain six simultaneous equation in six unknowns which 


are presented below in matrix form, 


idee O Q ILA, 8) 0 QO ay 2,0 
e,0l33603 1,0 -~2,0 ae 0 0 as -1,0 
BeILSO003 6,013363 qe -~2,0 Lgyay O) QO by @) 
~140759 2.918063 QO 1,0 -~2,0 1a b, 7 QO 
, 00019. ~ 140759 QO e) a0 -~2.0 by QO 
QO PeOOOM 1 Q e) 0 cae, by Q 


This matrix equation was solved by computer using a Gaussian 
elimination library routine, With the results obtained there- 


from, the coefficients of the process, K(z), were calculated 


as follows: Aig = 2 b, = .319530814 
So 2b) -li- bo = 1,70120542 
Dies j= 2), = by — b, = -, 39202846 
Cd 2b. _ b, ~ b, = -,59772472 
EN SS 2b, = b, = ~,03094094 
A(6) = -b) = -,00004212 


These process coefficients were supplied to the simulator 
program along with the coefficients of P(z) and Q(z) and 
the values of the four plant zeros to be included, The 
program then calculated D(z) for the single-rate sampled 
Gata system, For a test input of a noise-free ramp the 


output has no overshoot or inter-sample ripple and reaches 
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zero error in 0,4 seconds, See Figure 38, The response 
to a noise-free step is shown in Figure 39, where the 
"tuned" nature of the design is again exhibited, However, 
after the initial large overshoot the step response, too, 
is essentially ripple free and steady state in 0,4 seconds, 
For the noisy ramp signal shown in Figure 40 the sam- 
pled data system RE unstable, while for a step signal with 
the same amount of additive noise the response can be seen 
in Figure 41, Note that the time scale in Figure 41 is five 
times longer than in the noise-free experiments, and the re- 
sponse to the noisy signal shows no signs of settling down 
any further at the end of five seconds than it did at the 
@ma of the first second, Clearly, noise of this power cannot 
be tolerated in the input signal to a sampled data system 


regardless of the digital controller design, 
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8, Conclusicns, 

The behaviour of a unity feedback controi system has 
been simulate:i using a digital computer program, Using ‘nis 
Simulator orcorem studies have been made of digital contros 
methods for veraous plants, with oarticular attencion fFo- 
cused on the contrel of a Naval gun mount, The azinuth func- 
tion of the gun mount is described by a fifth order oi:fer- 
ential equation, and it is believed that the concivcions 
pertaining to it are applicable, at least in wart, to oth 
shipboard weapon systems, 

The use of multi-rate digital controllers was inves 
tigated, in every instance, the double-rate controller poro- 
vided faster response by one-half a sampling perivz, How- 
ever, only in tis special case of riople-free desiia For a 
step input was the double-rate response an unallioyen im- 
provement over the single-rate response, In gene: -., 
transient concroi and inter-sample ripple was not ictitee 
ably improved .y using a double-rate controller for the gua 


mount, For otixr plants of different frequency rfesvonse 
ae a, os a 


ay 


advantage mig.:t o@ taken of the fact that inter-samp.s 
ripple when using a Gouble-rate controller occurs «’* tiu21ce 
the error-samotang rate, 

Pulse tceansfer functions (or Z-transforms) 3:2 a funce 
tion Of the sawoling veriod, It was found that for certain. 
sampling rates valid simplifying assumptions coulis . aoe 
which facilitated the design procedure, Care snouiu oe 


+ 


exercised that these simplifying assumptions are mci. made on 


on 





the basis of a sampling rate which is too slow in relation 
to the damped natural frequency of the system concerned, 
In this regard the approximations made to the pulse trans-= 
fer function of the gun mount were inappropriate since they 
were based on a sampling rate roughly equal to the damped 
natural frequency instead of six to twelve times greater as 
must be the case for adequate control, 

Extreme sensitivity of sampled data systems to noise 
in the input signal was demonstrated, No satisfactory method 
of coping with input noise by design of the over-all process 
was found, So-called "variance reduction" processes showed 
little merit as the result of these simulation studies, The 
design criteria which came closest to producing satisfactory 
response to a noisy step input was the prototype with stale- 
ness factor of 0,5, 

Without compensation, either analog or digital, the 
gun mount is closed-loop unstable, Several digital control 
functions were developed which gave stable response to noise- 
free input signals. None of these gave stable operation when 
noise was present in a first order (ramp) input or a second 
order (acceleration) input, but for a noisy zero order (step) 
input digital control achieved a response that was at least 
bounded for the duration of the experiment, albeit unsatis- 
factory for most conceivable applications, It is conciuded 
that sampled data system sensitivity to additive noise in the 
input increases with the order of the input signal, 


The computer program described in Appendix If for 


17 





Cinding Z-transforms was very helpful in this investigation, 
It is hoped that this program, as well as the control system 
simulator upon the results of which this work was based, 

will prove useful to others in the field of sampled data 
systems, It is further hoped that this work has, in a small 
way, indicated the truthfulness of the assertion thnat digital 
compensation techniques introduce a degree of flexibility ana 
accuracy in a control system that is not presently achievable 


with state-of-the-art analog techniques, 


) 5. 





0, 
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APPENDIX I 


DIGITAL COMPUTER PROGRAM FOR SIMULATION OF A SAMPLED 
DATA SYSTEM 


General Description of Program ALSTAPP 
This program is a simulation of the hybrid control 


system shown below, 
DIGITAL. <—— ee Bune 
INFUT CONTROLLER, a2 OLD PLANT UT 
¥(t) q < 55% C(t) 


Z RAW ERROR D(z) PROCESSEO a H (s) Laie om 


E1(z2) ERROR E2(l2) yi Ce) 


1 
! 
J 
Cae aoe 












-—_——_—= = 


R(z) 
>. Ce. 
0 SS i erm (Ca peepee eee 

Three choices are available for the test input; namely, 
a unit step, a unit ramp, and a unit acceleration, The 
plant may be up to thirtieth order and is currently sub- 
ject to the following restrictions: 

(1) No plant zeros are allowed, i.e,, the numerator 
of H(s) must be a constant, 

(2) The plant may not have a pole at the origin higher 
than first order, i.e., limited to Type O and Type 1 servos, 

' (3) Only the constant term may be missing in the dif- 

ferential equation describing the plant, 

In essence, the plant is restricted to having a trans-= 


fer function of the general form 


VK 
Ge S$ (s+ AYS+ HR) (S* Ay) = _- Simple poles, Pi 


IK 
H(s)= — A(ket) sR + A(k) ok" + «-+ + AQ)S + A(1) 


or 


where A(l) is normally zero for Type 1 servo and is the only 


missing term allowed, 
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The rate at which the raw error, El, is sampled can 
be selected at will. This sample period T is chosen by 
the user and supplied to the computer, 

DELT is the small increment of time used in the simu=]= 
lation of the continuous plant by numerical solution of the 
plant's differential equation, This is done one thousand 
times, so that the total problem time is 1000 x DELT, 

DELT should be chosen as small as possible consistent with 
the desired total problem time, 

Although one thousand values of the solution are gen= 
erated and are available for use by a computer plotting 
routine, only every tenth solution is produced on the Ana-= 
lex printer, since it was felt that a column of one hundred 
values is sufficient for perusal by eye and more than ade= 
quate for plotting by hand, 

This program performs four major tasks, The first 
three are performed every time the program is executed, 
but it is optional with the user whether or not the fourth 
task will be done, These tasks ares 

(1) Simulation of a continuous, un-compensated, unity 
feedback system for the given plant, 

(2) Calculation of coefficients for a single-rate 
digital controller, 

(3) Simulation of the sampled-data system employing 
a single-rate digital controller using the coefficients 
calculated above, 


(4) Simulation, if desired, of a sampled-data system 
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employing a double-rate controller, The coefficients of 
the double-rate controller must be supplied by the user, 
and if they are not, the fourth task is not performed, 
Iinforination must be supplied to the program via data 
cards, Detailed instructions for the preparation and con- 
tents of these cards is given later, Seventeen (17) cards 
are required if the program is to do all four tasks, while 
only the first thirteen (13) need be supplied if simula- 
tion of the double rate controller is not desired, CAUTION: 
The cards must be presented to the computer in exactly the 
order and form specified, The slightest error will cause 


faulty operation or, perhaps, no operation at all. 


Discussion of Information Supplied by the User 

The general nature of the information that must he 
supplied to the computer is described next, 

(1) The user must choose the small time increment 
DELT to be used in the numerical method solution of the 
plant's differential equation, The smaller this increment 
the more accurate the solution, However, the total time 
simulated will be only one thousand times DELT, so the user 
should first determine the time span of interest...i.e., 
an estimate of how long it will take the system to achieve 
steady state, or how much of the transient condition it is 
desired to observe for lightly damped systems...then divide 
this total time by 1000 to obtain the value for D&LT, 

(2) The user must choose the sampling period to be 


employed in the sampled-data system with single-rate 
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controller simulation, This value is called T, 


(3) The user chooses the test input to be applied 
to the system, Unit values of step, ramp, and acceleration 
are available, 

(NOTE 1: These first three pieces of information, namely, 
DELT, T, and INPUT are supplied by data card number one, ) 

(4) The specifications of the continuous plant are 
supplied to the computer by data cards 2, 3, and 4, The 
order of the plant is denoted by K, This is the highest 
power of the Laplace variable, S, in the denominator of 
the plant transfer function, H(s). The numerator gain 
constant of H(s) is denoted by VK, The order of the pole 
at the origin is denoted by NTYPE and is limited to 0 or l, 
These values, K, VK, and NTYPE are supplied by data card 
number two, 

(5) H(s) must be known in the polynomial ferm, that 
is to say, a numerator gain constant, VK, over a denominator 
polynomial in S, The constant coefficients of this denomi- 
nator polynomial are denoted by A(i) and must be supplied 
to the computer, The only coefficient which may have a 
zero value is the constant term A(1). The number of co= 
efficients is denoted by NA and is supplied by data card 
number 3, Even though A(1) is normally zero (in a Type l 
servo) it must be included in the count, NA, on card 3 and 
its zero value shown in the first ten columns of data card 
number four (4), In the usual case of a Type l servo, the 


value of NA is one more than the order, K, of the plant, 
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(6) For the simulation of the sampled-data systems, 
a-transform techniques are implemented, Data reconstruc= 
tion is accomplished by a simulated zero order hold, which 
must be included in the Z-transform approach to the plant, 
If we call the plant and hold-circuit combination G(s), then 
G(s)= In E71 Gs) and the pulse transfer func- 

tion (or z-transform) is 
=r S z 
et): (ine) % AE. FS 
This pulse transfer function for the plant and hold must 


be supplied by the user, Its general form is a ratio of 


polynomials in i as shown below: 
=I =a —-NP 
eee Par Pare eRe Ue 
Q@) Q (4) ote Qa) z! + Q(3)27*+ coe OGg)iae St 


where NP is the number of terms (and also the highest 
power of g~*) in the numerator, P(Z), and NQ is the num- 
ber of terms in the denominator, Q(Z). In order for G(Z) 
to be a physically realizable pulse transfer function, the 
term Q(1) must be a constant, and in fact, is usually unity, 
(7) Although G(Z) is supplied to the computer in 
polynomial form, it must be Known in factored form, in 
order that any poles or zeros of G(Z) that lie on or out- 
side the unit circle in the Z-plane may be accounted for, 
Usually, the plant is open loop stable and only Zeros need 
be accounted for, but the program will also process infor= 
mation and provide the digital controller coefficierits 
When the plant has a pole outside the unit circle, All 


such zeros must be supplied by the user under the generic 
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term ZERO(i) and the number of them is denoted by NZ, 
Up to ninety=nine poles and/or zeros may be accounted for, 
(8) The overall pulse transfer function of the en-= 


tire system is denoted by K(Z) and is defined as 


«, SG) 
K (2) = Ris) 


For the purposes of this program, K(Z) must be a 
numerator polynomial only, This implies a "finite=-memory" 
process, If an infinite-memory or "non-finite settling 
time process" is desired, then it must be approximated 
by a finite number of terms of the infinite series ob= 
tained by long division, The program will operate on up 
to ninety-nine (99) such terms for K(Z),. The general form 
for K(Z) is 

K(z) = AKU) 2 + AR(2)B ++ + AK(MK)Z 
where NK is the number of terms provided, K(Z) must be 
supplied by the user, 
NOTE 2: In order to supply an array of numerical data to 
the computer two (2) data cards are required, The first 
card tells the computer how many values are in the array, 
and the next card (or cards) supplies these values, As-= 
suming for the moment that no array contains more than 
eight (8) values, the following information is supplied 


to the computer by the data cards listed below: 


Card Nos. Information 
3 & 4 Array of coefficients of plant 


differential equation (i.e,, 
denominator of H(s)). 


5 & 6 Array of coefficients of nua 
merator of G(Z). 
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Card Nos, Information 
7 & 8 Array of values of zeros and 
poles of G(Z) to be included 
in overall process pulse trans=]= 
fer function, K(Z), 


9 & 10 Array of coefficients of de= 
nominator of G(Z), 


lil & 12 Array Of coefficients of over=]= 
all process pulse transfer 
function, K(Z), 

(9) For certain investigations it may be desirable 
to limit the output of the digital controller to a signal 
of a given magnitude to avoid plant saturation, If so, 
the user must supply this value as a positive decimal 
value, E2MAX, and the program will limit the processed 
error from the controller to plus or minus this value, 

If no limiting is desired, a value of zero (0,0) should 
be given to £2MAX and the program then sets a magnitude 
limit of one million on the processed error, Data Card 
13 contains the data on E2MAX, 

NOTE 3: In very simple cases it may be that there are no 
zeros or poles of G(Z) that need to be included in K(Z). 
When this oceurs, obviously the datum NZ is zero and the 
array ZERO(i) does not exist, Under these conditions, 
prepare data card 7 as shown below 


Col. Ll Card 7 


and remove card 8 completely! 
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DATA CARDS for input to Program ALSTAPP 


Card No, Format Contents Remarks 


SS EEE 


DELT, T, INPUT Total/ time = 1000 

DELT.| Sample period 

T = V/ferror sampling 
rate, DELT is the ine 
ecrement used for simu=- 
lation of continuous 
plant, INPUT = 1] for 
a step, 2 for a ramp, 
3 for acceleration, 


K, VK, NTYPE | K & order of plant, 
VK = plant gain, 
NTYPE = type of servo 
(i j.ey, 1).2 239) 


é 
3 IT2 ! NA = no, terms in dee 
| | nominator of plant 
transfer function, 
Usually NA = K #1, 
H(S) = 
| VK. 
| | A(K+1)s* +. AC 2)S+#A(1) 


8F10.0 1A(I) IT=1, NA] A(I) is the array of 
coefficients of plant 
differential equation 
(i,e., denominator of 
H(S) above), 


NP = no, terms in nue 
merator of Z-transform 
of plant and zero-= 
: order hold whose Ze 
transtorm is 
! _ F(Z) 
AS ae =< 
P(I) I=1, NP P(1)z7"+P(2)2" 7400, 
! 1+0(2)27*40(3)2° "+... 
+Q(NQ)Z-NO#1L 
— 
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2F10,0,11 







I2,F18,0,I1 














. 


2 
8Fi0,0 
7 I2 NZ = no, zeros of G(Z) 
| to be included in over- 


all pulse transfer fca, 
K(Z). 


i 
ae 














ay, 





ri 





DATA CARDS for input to Program ALSTAPP (Continued) 


Card No. Format Contents Remarks 


8F10.0 ZERO(I) I=1,NZ = array of 
values of such zeros 
of G(Z) as are to be 
included in K(Z). 

I2 NQ = no, terms in dee 
nominator of G(Z), 

SF10,0 Q(I) I=1, No Q(T) = array of de= 
nominator cost ages. 


NK = no, terms in over- 
all pulse transfer 
function 
K(Z)=aK(1)27*+aK(2)277 
oo. TAK(NK) 25K 
SF10.0 

cess coefficients 

shown above, 

13 F10,0 E2MAX | B2MAX is the limit 
placed on the absolute 
value of the output of 
the digital controller, 


Notes to the User: 


AK(I) = array of pro-= 





The FORMAT conveys information concerning the mode and 
location on the IBM card for each item of numerical infor-= 
mation to be supplied to the program, If the following 
remarks are not sufficient, refer to the CDC Fortran Manual, 

"I" means an integer value, Whole numbers only, Do 
not show a decimal point, The integer value must be "right 
justified" in its allotted field, If the field is left 
blank, a value of zero will be inferred, 

"EP" means a decimal value, either less than cr greater 


than unity, A decimal point must be shown somewhere in the 
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allotted field, If the value is negative, the minis sign 
must also fall within the field, Precise location of the 
value (including sign, if any, and the decimal point al-= 
ways) within the field does not matter, 

The numbers following I and F give the number of col- 
umns on an IBM card that are reserved for the field of 
that particular datum, For example, take the FORMAT con- 
trol for data card number 2: Columns 1 thru 2 are reserved 
for the integer value of K, If the integer value of K is 
four, then card 2 should have 04 in its first two columns, 
The next eighteen columns (3 thru 20) on card number 2 
must contain 1500000,0 somewhere between columns 3 and 20 
inclusive, On both cards 1 and 2, column 21 must contain 
a single integer. 

On those cards having an I2 format, columns 1 and 2 
contain an integer value between one and ninety-oine, If 
the value is less than ten, it must appear in column 2 
(i.e,, the integer value must be right justified in its 
field, Decimal values need not be justified), In every 
case these I2 integer values tell the computer how many 
decimal values to read from the next card or cards, This 
brings us to the FORMAT 8F10.0, As before, the F1i9,9 
means a Gecimal value with ten card columns reserved for 
it. 8F10.0 means there may be up to eight such values on 
a single card (a card has eighty columns), and the com- 
puter will read another card if it has to, in order to 


read as many values as were specified by the previous I2 
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card, Most input arrays will contain eight or less values, 
however, and thus will fit nicely on a single cara, Note 
that these decimal fields of ten columns each begin in 
card columns 1, ll, 21, 31, ...etc, Card 1 has two decie 
mal fields of ten columns each, followed by an integer 
field of one (1) in column 21, Card 13 has a single 
decimal field occupying card columns 1 thru 16, 

If a decimal value of zero is to be entered, then 
place 0,0 in the appropriate field, For example, suppose 
array A has three decimal values A(1l) = zero, A(2) = 


twenty, A(3) = minus two anda half, Then Card 4 should 


contain 
Between Columns The Decimal Value 
lL and 10 OG 
iI and 20 20,0 
zl and 30 =2,3 


If no sign is shown, the value is taken as being 
positive, 

In order to use this program to test a system using 
a single-rate controller, only the first thirteen data 
cards are required, 

If you desire to test a system using a doublie-rate 
controller, then an additional four data cards ate re= 
quired, These cards contain information as to the double- 
rate controller constants and follow card 13 (no inter- 


vening blanks!) as shown below: 
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ADDITIONAL DATA CARDS FOR DOUBLE=RATE CONTROLLER 


SS eee eee eee 


Card No, Format Contents Remarks 


14 | 12 NAA WRA 12M the nas ener 
i in numerator of doublee 
rate controller, D(Z.). 
15 8F10,0 AA(I) I=1, NAA]| Array of numerator 
| coefficients, 
16 ie NBB 


NBB is the no, terms 
in denominator of 


é “abt +eet ae 


wee tAA(NAA) waned 


, —N + 
oo» BB(NBB)Z oe 


17 | 8F10,0 BB(I) I=l, NBB] Array of denominator 
| coefficients, 


Example Problem using Program ALSTAPP=: 








H(S) = ag SOO 
S~ + 955° + 1750S 
By looking at H(S) we can see that VK = BO000_90 
ae 
NA = 4 
A(1l) = 0,0 


A(2) = 1750.0 
A(3) = 95,0 
A(4) = 1,0 
NTYPE = 1 


Let us choose to observe the response of this pliant for 


one second after the application of a unit step input, 


ai 


So DELT 1/1000 = 0,001 


INPUT l 


fl 








Let us select a sampling rate of 160 times per sece 
ond for the sampled data system, i,e., T= 9,01, 


We prepare the first four data cards as shown belo 


Column Column Column Column 
2 eee 21 a 
Card 4 / j0,0 »4750,0 95.0 + O ais 4, 
ee Ne ——- 
Cara 3) 304 | | 7 i Card 3 
ce rer anor oes es, 
Card 2h j03 80900,0 11 Card 2 
CS ee 
Card lp’ ,O,001 yO, O1 iL \ Card 1 
i 
i J 
By partial fraction expansion and tables of Z-trans= 


forms we find the pulse transfer function of the vLlante 
and=-hold to 5# 


in factored form: 


G(Z) = 01062 = + Oy 20932727") 1 + 2 296727527) 


(1-2 = a eee a 1) (a - 0,4965852" 5} 


coz) = 20196007" + ee 033672277 + ,0065842" ae peepee: MY, 
aes go) + 1.662127Z°2 = 0,3867412°° 


There is ons zero outside the unit circle at -2,367275, 
Choesing a minimal prototype process we Find that 


=~] =2 


K(23) = 0, 2520022 + 00,7479 382 


and we seé the tollowing data: 


NP = 903 P(1) = ,010600 
P(2) = .033672 
P(3) = ,006584 
NZ = O1 ZERO (1) = =2,967275 
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NQ = 04 Q(1) 
Q(2) 
Q(3) = 


Q(4) = 


AK(1) 
AK( 2) 


NK O2 


So we prepare data cards as fo 


Column Column 

1 11 
Card 12/10. 252062 0.747938 
Card 1l 02 | 
Cara 10% j1.0 1-2, 275386 
Card 9 / ,04 | 
Card 8 =2,9667275 | 
cara 7 |“'01 | 
Card 6 “1 010600 | 033672 
Card 5 103 





120 
—-2.275386 
PeCoz 27 
-0, 386741 


0.252062 
0,747938 


llows: 


Column 
21 


1) 662127 


1+,006584 


Column 


\-0, 386741 


| 
oe 
| 
l 










Suppose we do not desire to limit the output of the 


digital controller, 
card 13 as shown below: 


Column 
1 


va 


(0,0 


If it is not desired to simulate a double-rate digital 


Then set E2MAX = zero by preparing 


Card 


controller, then only the first thirteen (13) data cards 


are required, 


If double=-rate action is desired, then the user must 
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supply the coefficients of the double-rate controller, 
The pulse transfer function of the double-rate controller 
has the following general form: 

AA(1) + AA(2)Z5> + AA( 3)Z5° + ... + AA(NAA)Z 


D(z.) = i i 2 
BB(1) + BB(2)Z5- + BB(3)Z5° +... + BB(NBB)Z,\>P** 


———_ 








where all coefficients have been normalized with respect 
to BB(1) and BB(1) = 1,0, 

The format for cards 14 thru 17 is identical to that 
for cards 5 thru 12, 

In the double-rate controller the error-sampling 
period, T, remains the same as it was for the single-rate 
case, It is the output from the controller that occurs 


at the double-=-rate, 


Modification to allow for addition of white noise to the 


ai TT St Ca = <eeeceeCCOO S 


test input, 
When used as described in the main body of this 


appendix, the program simulates a control system whose 
input is completely deterministic, Gaussian noise of 
selected mean and variance may be added to the input by 
making an addition to the first data card and inserting 
a new data card between the previous card 1 and 2, as 
follows, 

1. The presence of noise in the input must be ree 
quested by placing the integer 1 in column 31 of the first 
data card (which already contains DELT, T, and INPUT). 

2. The new data card, which must follow card 1 


above, contains an integer number between 01 and 99 right 
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justified in columns 1 and 2, and a decimal value in 
columns 3 through 20, The decimal value specifies the 
selected mean value of the noise, and in most cases will 
be shown as 0,0, 

The integer number specifies the number of uniform- 
ly distributed random numbers from which each normally 
distributed random number is obtained, For a good approx- 
imation to white noise, this number should be between ten 
and twenty, The variance of the Gaussian distribution 


depends on this number, N, as shown below: 


1 lL 
N° 12 


Variance = 

As N approaches one, the noise approaches a uniform 

Gistribution between =-1/2 and +1/2, for a selected mean of 
zero, 


The remainder of the data cards are then included 


in exactly the order and format described earlier, 


Modification to allow the program to use single-rate digi- 
tal controller coefficients supplied by the user, 

When used as described in the main body of this 
appendix, the program computes the coefficients of the 
single-rate controller using the following formula: 

D(Z) = atgy » pS, 

In order to do this, the user must supply data on 
the numerator and denominator of G(Z), the overall pulse 
transfer function, K(Z), and which zeros, if any, of G(Z) 


that are included in K(Z), CAUTION: Slide=-rule accuracy 
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is usually not sufficient except in very simple cases! 


It is recommended that a desk calculator be used to pro-= 


vide this data accurately to at least five decimal places, 


In case the user wishes to calculate his own cone= 


troller coefficients for the single=-rate sampled-data 


system simulation, he may cause the program to omit the 


calculation of D(Z) and to employ the controller supplied 


by the user, In order to do this, the data cards are ar= 


ranged as shown below: 


Card No, Contents 

Ht DELT, T, INPUT 

2 K, VK, NTYPE 

3 NA 

4 A(I) I=1, NA 

5 Blank card! 

6 NAA 

7 AA(I) I = 1, NAA 
8 NBB 

9 BB(I) I = 1, NBB 
10 E2MAX 
11i* NAA 
12* AA(I) I = 1, NAA 
13* NBB 


Remarks 
Same as before 
Same as before 
Same as before 


Same as before 


Indicates D(Z) supplied by 


user 


Single-rate controller 
efficients 


Single=rate controller 
efficients 


Single=rate controller 
efficients 


Single=rate controller 
efficients 


Same as before 


Double=rate controller 
efficients 


Double=-rate controller 
efficients 


Double=rate controller 
efficients 


dil 


co= 


co= 


co= 


co= 


CoO= 


co= 





Card No, Contents ReEMacks 


14* BB(I) I= 1, NBS Double=rate controller coq] 
efficients 


*Optional, 

The first four cards are just the same as they would 
be for a deterministic input where D(Z) is to be computed 
by the program, A blank fifth card signals the program 
that D(Z) will be supplied by the next feur cards, These 
next four cards (6, 7, 8, and 9) give the Gata for tne 
single-rate controller exactly as described eariier in the 
main body of this memo for the double-rate controiler co= 
efficients, Card 10 gives SzMAK, If double=-rate action 
is also Gesired, the last four cards must give the couble=- 


rate D(z). 
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PLOW DIAGRAM for Frogram ALSTAPP 
START 


Dimension all arrays 


Read DELT, T, INPUT, NOISE 


Branch on NOISE NOISE # 0 


NOISE=0 ‘ 


\ 


Read K, VK, NTYPE 


Read NA, (A(I)I=1,NA) 


Read NRAN, BIAS 











Compute theoretical 
standard deviation 





Read NP 


Branch on NP 
VE 6 





Read P(I) I=1,NP 


Generate NRAN ran= 
Read NZ 


Gom numbers of uni- 
form distribution 
between O and l 
(ant on 
Read ZERO(I) 
NZ=0 


Form a Gaussian 
random variable, 
RANO(J) from them 


J=1000 


J<1L000 


Read NQ,(Q(I)I=1,NQ) 
Read NK, (AK(1I)I=1,NK) 


| 


Read E2MAX 









Compute mean of 
generated noise 






aM 









Print heading, plant type, 
gain, and differential equation 


Print D(Z) as Print K(Z) 
supplied by user 


Print P(Z) 
(mes 
Print Q(Z) 
Compute MSWITCH 
Print ZERO 
array 


Compute RATE 
Print RATE 


Branch on B2MAX E2MAX=0 


Set E2MAX= 
one million 









Print E2MAX 





i 


Print type of INPUT 


NOISE = 0 Branch on Norse )—NOISE # 0 


Clear the array of Print standard 

RANO(I) values deviation and 
measured mean 
of input noise 


Print heading for results of 
continuous, un=-compensated, run 


G9 








Set all states of the 
plant equal to Zero 






i 


Time (1)>= 0 






INPUT = l INPUD = 2.00 


R(1)=0,0+RANO(1) | 


- Branch on JJ ee Gas 
JJ = 2 





R(1) = 1.0 + RANO(1) 


Branch on INPUT 





so 





G59 


Compute plant response 


by solving differential equation 
for forcing function U 


Time (N+1) = Time (N) + DELT 


Branch on INPUT NE SE 


INPUT = 2 


INPUT = 1 


R(N+1)=1+RANO(N+1 ) R(N+1 )=TIME(N+1) +RANO(N#41 ) 


R(N+1 )=TIME(N+1) “+RANO(N#1) 
\ S 

C40) : Yes Is N less than 100? 

No 

Is N less than 1000? ac 


Branch on JJ 3 


Yes 














Call GRAPH rou] 


tine if avail= 
able 


JJ=1 





Print results 
of run 


Branch on JJ 


JJ=3 
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ee 


Branch on NP NP #0 












Compute single rate 
controller coeffic= 
ients, given G(Z), 
K(Z), and the list of 
zeros of G(Z) to be 
included in K(Z), from 


ee (072 
DZ “Seeay trata 


Print computed filter 
coefficients 






Initialize the 
digital controller 


e- 


set M = MSWITCH 


Lf 


Compute PARTSUM = PD (AA (I)* EIN(I) - BB(I) * EOUT(TI)) 





Sense E1(N) = R(N) —- Y(N) 


Compute EOUT (1) = AA(1)* EIN(1) + PARTSUM 


1 


= Is |HOUT(1)| less than E2MAx? )—*eS 
Set B2(N) = B2MAX Set E2(N) = 
with proper sig BOUT (1. 


co) ~ 


' 


U = VK * E2(N) 
Shift filter number sequences 


LL.7 












Read NAA 





Branch on NAA 









Read AA(I) I = 1, NAA 
Read NBB, (BB(I) I= 1, NBB) 


Set remaining AA(I) and BB(I) to zero 


= lor 2 Branch on JJ 


JJ = 3 


Print double-rate 
controller coefficients 


Compute M2 


uw 
ma 


Set M = M2 
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NEO = 1 Branch on NEO Seo 


Compute PART1 = D> AA(21I)* E1N(TI) EL(N) = R(N) = Y(N) 










Set E1(N) = 0,0 Compute PART] = 


2 A(2I-1) * EIN(I) 





Shift input sequence 





Set NEO 


Set NEO = 2 
* a 


Compute PART2 = 7 BB(I) EOUT(I) 


EOUT(1) = PART1 = PART2 


No 


Yes 


Is |zouT(1)| < =2mMax? 


Set 52(N) = s#uT (2)| 
) 


U = VK * E2(N) 


Shift output sequence 
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APPENDIX II 


DIGITAL COMPUTER PROGRAM FOR FINDING Z=-TRANSFORMS AND 
ROOTS OF POLYNOMIALS 


General Description of Program Stappl 


This program is a combination rootfinder and Z=-trans= 
form taker, It will operate in one of three modes, as 
specified by the user, These modes are: 

Mode 1 

Program operates as a rootfinder only, The user sup= 
plies the degree and coefficients of the polynomial whose 
roots are to be found, Handles polynomials up to and in-= 
cluding thirtieth degree, 

Mode 2 

It is desired to take the Z=-transform of a transfer 

function which is not known in factored form, such as 


H(s) = VK 
a(N+1)s‘t+a(n)sS-t+ . |. + A(2)s + ACL) 


The program first finds the roots of the denominator 
polynomial, Then it takes the Z-transform of the above 


plant transfer function preceded by a Zero-order hold, i.e. 


-st _ =] 
G(Z) = x, (les oH(s)} = (1-2 ) H(s) a P(Z) 


and presents the z-plane poles of G(Z) as well as the 
numerator and denominator polynomials in gt. Then it 
operates as a rootfinder on the numerator, P(Z), in order 
to display the zeros of G(Z). 

Mode 3 


It is desired to find the Z=-transform of a plant and 
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hold combination when the transfer function of the plant 


is known in factored form, such as 


VK 
s S-P) S-P. S-P, crea S=P. 
As in Mode 2, the Z-transform is taken with the plant 


preceded by a zero-order hold, i.e, 


ez) al-z 4) Als) . B(z) 


The program displays the Z-plane poles of G(Z) and the 
numerator and denominator polynomials, It then operates as 
a rootfinder on the numerator polynomial, P(Z), in order to 


display the zeros of G(Z), 


Restrictions, 

Mode 1 is relatively unrestricted, It will find the 
roots (real and imaginary parts) of polynomials in the 
following form (for N up to and including 30); 

P(x) = A(N+L) x™ + a(x" 7+ + 22. 4a(2)x + ACL) 

Modes 2 and 3 are restricted to polynomials of degree 
ten or less, with all roots lying on the real axis and no 
repeated roots, Mode 2 will accept a polynomial whose 
eleventh root is at the origin, In both modes the Z-= 
transform operation is predicated upon the plant being a 
type 1 servo (i.e, the Laplace transfer function has a 
single pole at the origin), 

Since the program will not take the Z-transforin 
when roots are complex, Mode 2 checks the magnitude of 
each imaginary part with respect to the corresponding real 


part of each root it finds, If the magnitude of any 
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imaginary part exceeds 1/100th of the corresponding real 
part, the program stops, If all imaginary parts of the 
roots found in the first phase of Mode 2 are less than 
1/100th the magnitude of their corresponding real parts, 
then they are considered negligible and the root is treat= 
ed as entirely real, The program then proceeds as de= 
scribed earlier, 

If, thru error, Mode 2 is selected when the unfactored 
denominator of H(s) is higher than degree ten (exclusive 
of root at origin) then the program stops after finding 


the roots of the denominator, 


Mathematical Methods: Root-finder 

The rootfinding techniques is a modification of Bair- 
stow's iterative procedure for finding a quadratic divisor 
of the form x°—PP x-QQ. When this divisor is found suf-= 
ficiently closely, it is solved by the quadratic formula, 
and divided out to leave a reduced polynomial two degrees 
lower, The criterion of suitable closeness is applied 
first to the change, DPP, in PP and then, if PP is suf-= 
ficiently accurate, to the change, DOQQ, in QQ, In both 
cases it is considered sufficiently accurate when the 
magnitude of the next correction to the coefficient is 


10 of the magnitude of the coefficient itself, 


less than 10° 
The essence of Bairstow's iterative procedure is to 

divide a polynomial, say A(x), by the trial quadratic 

divisor to obtain a first quotient polynomial, say Q(x), 


olus a first remainder, Then this first quotient is 
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Givided again by the same trial divisor to obtain a sec= 
ond quotient, say T(x), plus a second remainder, Then 

the corrections necessary in the coefficients of the 

trial quadratic divisor are determined from a suitable 
combination of the coefficients of the first and second 
remainders, A symbolic example will illustrate the method, 


Let F(x) = A(N+1)o“ea(n)o t+ 5. 4a(2)x#a(1) = the poly- 


nomial whose roots are to be found, 
Let the trial quadratic divisor be TD(x) = x7=PP, x-00 


Now F(x) = TD(x)*[Q(N+1x8~ +0(N) x87 A+... +0(4)x+0( 3) | 


+Q(2)xtQ(1) 


where the first quotient is Q(x) inside the brackets and 
the first remainder is Q(2)x+Q(1l). 


Dividing again by the trial quadratic divisor: 


Q(x) = TD(x)>['T(N-1) x87 F47( 9-2) 24+, +704) x47( 3) | 
+T(2)xt+T(1) 


where the second remainder is seen to be T(2)x+T(1) 


The corrections to the coefficients of TD(Z) are computed 
from the following formulae: 


M = PP,T(1) + QQ.T(2) 
= T(1)2-M, (2) 


mene t(2)-O(1) 7 T(1)-Q(2) 
pog = Me2(2) = 7(1)-0(2) 


The modification of Bairstow'’s procedure came about as 
the result of several runs using test polynomials, It 
was found that convergence in the case of non-repeated 
roots was more consistently obtained if each set of iter] 
ations on successive reduced polynomials was begun from 


the origin, In case there are repeated roots the return 
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to the origin is delayed until all roots at that location 
have been found, It was also found that an occasional 
set of iterations would not terminate even though no fur- 
ther refinements were being made in the coefficients of 
the trial divisor, When this occurs the program assumes 
the trial divisor is satisfactory after fifty iterations, 
No degradation of accuracy has been observed due to this 
artificial convergence device, 

Because of the sharpness of the criterion for a 
satisfactory trial divisor it appears that accuracy of 
root values to at least seven significant figures can be 
expected, Since the sum of the roots of an Nth degree 
polynomial is the coefficient of the term containing 
the variable to the (N-1)th power, the program calcu-] 
lates a check sum from the root values found and dis- 
plays this check sum along with the true value as an in= 
Gication of the accuracy achieved, This also provides 
the user with the means to tell at a glance whether or 


not the procedure converged, 


Mathematical Methods: Z-transformer 

The Z=-transform portion of the program is restricted 
to type 1 servos preceded by a Zero order hold, That is, 
the general form of the function whose Z-transform is to 


be taken must be 


-sT 


l=-e s 


VK 
etsy = se S(s-P,)(s=-P.)...(S=P_) 
Ss s(s Pi Ss Po eee SS Py 


Where s is the Laplace variable, T is the sampling period, 
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and N is the number of poles on the real axis; not in- 
cluding the pole at the origin, 

Note: 

All poles must be real and none can be repeated, 


By definition, z=esl so by partial fraction expansion we 


can write 
G(Z) = th, BAS Bs C(1) oy C(N) 
s-P, ~ S=-P °°° s-Pn 


The program first computes the partial fraction co-= 
efficients, A, B, and C(i) i=1, 2, ..., Ne The Z- 


transform of each term is then taken as follows: 
A ATZ 
i wer 
Ss (Z-1) 
eer 
Ss 2-1 


Se ) etal) 


S - Pi ae arit 





The poles of G(Z), exclusive of the one at unity, are 
seen to be at Z= eri for i = 1; 2;,°23.,0" 

The program then proceeds to combine the individual 
Z2-transforms over a common denominator, cancel the external 
-1, 


factor (1-2 , and multiply through by the plant gain con- 


Stant, VK, 
The final results are displayed as the coefficients 


of the numerator and denominator polynomials shown belows 


G(z) = P(1)z + p(2)zN-1 4+ p(3)gN-? +... + P(N)Z + P(N41) 


o(1)28F* + e¢2yz™ + a¢apz%h* + 0. + a(ntl)z + O(ne2) 


or, equivalently: 
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G(Z) = 





Discussion of the Data Provided by the User 

The general nature of the input data for each mode 
is discussed below, The sequence in which the data are 
presented is vital, and is correctly shown in the discussion. 
Mode 1: 

1) MODE = the integer number l, 

2) N = an integer value between 1 and 30 representing 
the degree of the polynomial whose roots are to be found, 

3) A(i) = a linear array of decimal values repre=-= 
senting the coefficients of the polynomial whose roots are 
to be found, These values must appear on the data card(s) 
from left to right in the same order as they are in the 
polynomial when written in descending powers of the vari= 
able, Example for N = 4; 

a(5)x* or a(4)x°> a5 a(3)x* + A(2)x + A(1) = F(x) 

Mode 2: 

1) MODE = the integer number 2, 

2) N = an integer value between 1 and 10 representing 
the degree of the denominator polynomial of the plant 
transfer function, H(s). 

3) A(i) = a linear array of decimal values repre- 
senting the coefficients of the denominator of H(s) in 


descending powers of A, Example for N = 5:3 





H(s) = 5 


& 2 


a(6)s> + aA(5)s° + A(4)s° + A(3)s* + A(2)s + ACI) 


135 









' 
-_ aay @ Ps | 
afte Die ———ow fom et > 
— A TT | ea | om ee om 
oo za = © & ae 6 =—& 
—_— — ee ieee —_ 
= «= & 6 =D; ee oe ee eae 
+ —_— —_ A « we 


2 





For a type 1 servo, A(1) = 0.0 and N is allowed to 
be as high as ll, The pole at the origin will be removed 
by the program and N will be reduced by 1 when a Zero 
value is read for A(1l). It is optional with the user 
whether or not he removes the pole at the origin himself, 
or allows the program to do it, 

4) PERIOD and VK = both are decimal values, PERIOD 
is the reciprocal of the sampling rate for which the Z-= 
transfer is to be taken, VK is the numerator gain con= 
stant of the plant transfer function, H(s). 

Mode 3: 

1) MODE = the integer number 3, 

2) N = the integer number of poles of the plant 
transfer function, exclusive of the pole at the origin, 

3) P(i) = a linear array of the decimal values of 
the poles of H(s). Does not include the pole at the origin. 
All poles must be simple and must lie on the real axis of 
the s-plane, 


4) PERIOD and VK = same as in Mode 2 discussed above, 


Detailed Instructions for Preparing Input Data 
All input data are provided to the program through 
the use of IBM cards, The "normal" number of data cards 


required by each mode of operation is tabulated below: 


Mode No, of data cards required 
L 3 
Z e 
3 4 
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"Normally" it is expected that the number of decimal 
values in the linear array will be eight (8) or less, so 
that the entire array can be contained on a single card, 
For arrays of more than eight values, additional cards 
are used as necessary, 

For convenience, the following table summarizes the 


requirements for input data: 


TABLE A 











— 
: 


x 
x 


* Use additional cards 3 as necessary for arrays contain= 






= 





ing more than eight (8) values, 

The "Format" information needs to be explained, The 
format letter I means that the value is an integer, [In= 
tegers must be shown without any decimal point and must 
be “right justified" in their allotted space, Right 
justified means that the right-most digit of the integer 
value must appear in the right-most column of the space 
allotted for that value on the card, This applies only 
to integer values, 

The format letter F indicates a decimal value anda 
decimal point must be used, If the value is negative, 


the minus sign must appear, Unsigned values are taken 
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to be positive, The decimal value, including the deci= 
mal point always and minus sign if necessary, may appear 
anywhere within the allotted space (or "field"), 

The number following the letter I or F is the number 
of columns on the IBM card forming the allotted space or 
field for that value, There are eighty (80) columns on 
a card and they are assigned from left to right in order, 

When a number precedes the format letter, it tells 
how many times that particular field is repeated, Thus 
card 3 contains eight identical fields of ten columns 
each, The next card may be read by the program as a con= 
tinuation of card 3 if more than eight values are in the 
array. 

Card 1 uses only the first column, This column must 
contain either the digit 1, 2, or 3. 

Card 2 uses only the first and second columns, If 
the value of N is nine or less, then skip column one and 
place the digit in column two, 

Card 3 uses as many identical fields of ten columns 
each as are necessary, Note that these fields begin in 
columns 1, 1l, 21, 31, etc. 

Card 4 uses the first ten columns to contain the deci- 
mal value for the PERIOD, The next twenty columns (i.e., 
11 through 30) are reserved for the decimal value of VK, 
Decimal values may lie anywhere within their allotted 
field as long as the decimal point is shown in the field 


and the minus sign is in the field for negative values, 


138 





Let us now show the preparation of data cards for 
some representative problems, 
EXAMPLE PROBLEM - Mode 1 | 
Find all the roots of 
F(x) = x/+ 16x°+ 102x°+ 371x4+ 838x°+ 1200x°+ 1012x+ 420 
Since all we need is a rootfinding operation, we choose 
Mode 1l, 
The degree, N, of the polynomial is 7, The polyno-= 
mial is written in descending powers of x so the three 
data cards would be 


Column Column Column Column Column Column Column Column 
1 11 21 31 41 51 61 71 


ja 16,0 102.0 371.0 838.0 1200.0 1012,0 420,0 (Card 


2 
Il a. Card 
| z 
The answers are: 
Roots 1 and 2 = -1,.0 + j 1,0 
3 and 4 = -1,5 + j 1.658312 
5 and 6 = -2,0 + j 1.414214 
7 = -7,0 
EXAMPLE PROBLEM - Mode 2 
Find the Z-transform of the following configuration 
Zero-order 
hold plant 
1 -sT eee 4 
—__———_j—jEs= Hs) — 
20 samples/sec 20 samples/sec 


where the plant transfer function is 
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1,6064235 x Or 


s” + 134s* + 5064s° + 732745 


H(s) = 5 
+ 356983s 


Since H(s) is not known in factored form we choose 
Mode 2, 
The degree, N, of the denominator of H(s) is 5*, 
The denominator of H(s) is written in descending 
powers, Note that A(1) is missing! 
The numerator gain constant, VK, is 1.6064235 x 10° 
The sampling rate is twenty samples per second so 
the sampling PERIOD is .05 seconds, 


The data cards are prepared as shown below: 


Column Column Column Column Column Column Column Column 


1 Ly 21 Se 41 SL 6l 71 
0 mame aT Tp a 
We OS 1606423500,0 
i 
l1.0 134.0 5064,0 73274,0 356983.0 0,0 








* User may remove pole at origin, if desired, and enter 
4th degree denominator where A(1) = 356983, 
The answers are: 


Z-transform = 


- P(5)Z7> 
=4 =5 
ea) 2. +. 0(6)2 
where 
P(l) = 1,54749758 OO) = 0 
P(2) = 16,74877866 Q(2) = =2, 33676594 
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oo] =D a3 a4 
! + - Z ~+P(4)Z “+ 
o(z) = PUZ) = PUL)2 +P(2)27*+P(3) P(4)z 


Q(Z) Q(1)+0(2)271+0Q(3)27740(4)Z7 3+ 





P(3) 
P(4) 
P(5) 


16, 29465286 Q(3) = 1,92220653 
2.04809669 Q(4) = = 67242310 
- 901866472 Q(5) = - 08821342 
Q(6) = = ,00123091 


The zeros of G(Z) are at Z = = 00988230 
~ ,13606251 
- ,91925736 
-9,75793506 


The poles of G(Z) are at Z 0,57694981 
0,427414983 
0, 31663673 


0,01576442 


EXAMPLE PROBLEM - Mode 3 


Find the Z-transform of the following configuration: 


Zero—-order 
hold olant 


=} 





a 


20 samples/sec 20 samples/sec 


where the plant transfer function is 


ity = 1,6064235 x 10° 
~ gs(s + 11) (s + 17 s + 23 s + 83 


Since H(s) is known in factored form we choose Mode 3. 


The number, N, of real poles (not counting the one 


the origin) is 4, 


9 


at 


The numerator gain constant, VK, is 1,6064235 x 10°, 


The sample PERIOD is .05 seconds, 
The data cards are prepared as shown below: 


Column Column Column Column Column Column Column Column 


1 11 21 31 41 51 61 far 
—— 
i205 1606423500,0 
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Column Column Column Column Column Column Column Column 





1 1l pal 31 41 51 61 71 

 —— — ——————EE——— 
fee =17.0 -23.0 -83.0 Card 
3 
lo4 Card 
2 
i i 

\3 Card 
a ee eee 


The answers are the same as those obtained in the example 


problem for Mode 2, 
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FLOW DIAGRAM FOR PROGRAM STAPP1 
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Read MODE 
Mode = 1 Branch on MODE Moe: = 
Y Mode|j= 2 ’ 
Read N 
Read PTREAL 
(i) for™ ie= 
Y l1 to N 
Read A(K + 1 —- I) Read A (K + 1 = I) 
Per I = 1 to K for I = 1 to K : 
()- 
Yes Is A (1) = 0? 4 
Call Call 
Roots (N) No ZTRAN (N) 
Remove the pole 
a at the origin ! 
Stop A (K+t1l-1I) = P(T) 


for I= 1 to N 


* Call 
Roots (N) 
Call 


Roots (N) 


No 
No Do all roots have Yes 
negligible read parts? 
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APPENDIX III 


Determination of stability of a third-order 
System as sampling rate is decreased 








G2) = K(2 +2+5_,2 where K= £0,000 
(F JE ($s . . a As 57./42°x10 





i” Be ~3.7020x/0°5 
Q(z) ou (2-1)4 t) = 2z-e a= 6 D= ~0.48535%105 
Multiplying out: 

3 or | get zest ar =957 -FST 
Qia)= z-[é"+e +1 ]2 + [e +6 \ere 


Pla) = K[AT-6(€0 +e +1)- clare”) - D(are"") | 2 
+ K[-at(e*5 2°") + B/E +e" en") + C(1+ 22°") + D(te2d ies 


+ ee Bens Ge pe 





~—_________—- © (s) 






C(t) 





80,000 
S(stas \s+70) 






3 2 3 
Characteristic equation: A(2)= P(2)+ Q(2)= 432 + Q,2+Q,2+Q,y 
where Qa=1 
= Ya) —AST _ 
Qo = KATA KG EO KBe” 2KC~KCE ~ 2KD- KDE — KB-B SPT 


-A2ST 70 = = = an JOT: 
ne age — KATe’ + Kee" "+ kKGe”'. Kye ae } ee 
+KD + 2KD eae eit re San e IST 
-9ST ~9ST -70T are 
Mo = KATe -KBe —~KCe - iS - € 


2 





Note that G(z) has been expressed in positive powers of 
z since we anticipate making a change of variables, Later 
on, z-transforms are expressed using inverse powers of z 


in order to correspond to the concept of difference equa-= 








tions. 
. L+w 
Using the Bilinear Transformation: Fas Peer 
3 a 
= |/+wW ated iatend 
A(z) = As (1 *) + Ag (2) + oe (#2) + Lo 


= Oe [4 + Sas + 3a? | 4 Az | I cam ay co | J Aa, |1-40~-c"+ ws [+ Ga [1- Beat 30% os 


_ tw [s-Ggtd,-Ge| i ts" [3G-Ga- 4, + Fie] + ts] 33 + 42-4, 34, | + [4s+4,+4,+%] 


Routh=-Hurwitz criterion for third order equations 
all roots in left half plane if b, b, = b, b, pO. 
where bo = G3 + Ga+ Gy, + Go 
tee Sa, + Q,-2,- 34, 
bn = 34,-42-4, +3 
bo = Az ~-42 +4, -4o 
In terms of the system parameters, 
bee SS VKAT) +e" (KAT) +S KAT) + (KAT ~ KE ~ KP) 
by = © (3KAT+2KG+2) + e” "(Kat -2KG-2) 
pS" (KaT-2K8-2) + (KAT —-3KC- 3KD +2) 
b, = oe (3kat -4kB- 4) + e”"( KAT - 4«C) 
+ Ee (KAT-4KD) + (KAT + KO +KD + 4) 
e@ *T(KatT) re (-KAT + 2KB + QKC + 2) 
+ ET (eat + QKB tT AKD +2) + (KAT + SKC + 3KO +4) 


b3 


A computer program to solve the above equations ose 
one hundred values of T and the results obtained are shown 


on the following pages, 
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wewOB STAPP2 
PROGRAM STAPP2 


o DETSRMINATION OF ROUTH=HURWITZ CRITERION FOR THIRD 
ORDER SYSTEM 

< AS SAMPLING PERIOD IS VARIED 
VK=80000, 


A= .000571428 
B==., 000031020 
C= .000035555 
D=—,000004535 
T=0.0 
N=0 
PRINT 199 
199 FORMAT( 32H PERIOD ROUTH CRITERION) 
100 T=T+,001 
BXP25=2,71828**(~25,0*T) 
EXP70=2,71828**(=-70,0*T) 
EXP95=2,71828**(-95.,0*T) 
AO=EXP9 5*( VK*A*T=VK*Be-1l, ) —EXP70* VK* C=EXP25*VK*D 
OAL=EXP95*(VK*Bt1, )+EXP70* ( —VK*A*T+VK*Bt2, *VK*¥Ct1, )+ 
1 EXP25*( =VK*A*T+VK*Bt2, O*VK*Dt1, ) +VK*C+VK+D 
OA2=EXP70*( =VK*B-VK*C-1, )+EXP25*(=VK*B=-VK*D-=-1, ) +VK* 
(A*T=2, *C=2, *D) 
1 =] ,0-VK*B 
A3=1.0 
BO=A3+A2+A1=A0 
B1=3,0*A3+A2=Al=-3,0*A0 
B2=3,0*A3—A2-Al1+3,0*A0 
B3=A3=A2+A1=A0 
ROUTH= B1*B2=BO*B3 
PRINT 200,T, ROUTH 
200 FORMAT(F10,4, F20.5) 
N=N+1 
IF(N=100)100, 300, 300 
300 STOP 
END 
END 


(FORTRAN program for the CDC 1604) 
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PERIOD ROUTH CRITERION PERIOD ROUTH CRITZRION 


~0010 ~ 900060 ~O0510 =-2, 33345 
~ 9020 .00418 els 210. ~2,44305 
~ 9030 ~01226 ~Us30 =2, 29203 
~ 0040 202524 ~ C540 =-2,66273 
~ 9050 ~ 904279 »9550 =2./7266 
~0060 06411 ©9560 =2,88256 
~9070 - 08814 ~9570 =2,99239 
, 0080 Pe be or ~ 9580 -3,10199 
~ 0090 2 a760 ~0590 =3,21141 
~ 0100 - 16483 - 9600 =-3,32057 
~90110 - 18818 ~ 0610 =3,42940 
.0120 e 20579 ~ 9620 =-3.05/57 
~ 0130 JeZool ~ 9630 -3,64593 
~0140 - 23874 - 0640 =-3,75354 
~O150 - 24688 ~ 90650 =3,86066 
~9160 - 24986 ,- 0660 -3,96725 
3070 © 24736 2067/0 =4 07329 
eOilisl®, Je agio - 0680 =4.,17874 
~9190 wee oee ~ 90690 =4, 28358 
~9200 - 20541 ~9700 =-4.,38776 
~0210 ~17977 ©0710 =4.,49128 
~9220 ~ 14839 20720 =4, 59410 
~Olzs0 JAS 7 ~0730 =4,69620 
~9240 . 06889 ~90740 =4.,79756 
—» ,0250 STABILITY,02112 ~0750 =4,89816 
~ 9260 LIMIT =—-,03171 ~0760 =4.99798 
~0270 =, 08939 ~90770 -5, 09700 
~0280 ~,15167 ~9780 =-5.19522 
~0290 =, 21831 ~0790 =9, 29261 
~ 0300 =, 28905 - 0800 =5,. 38915 
~9310 ~. 36365 - 0810 =-5.48485 
yOs20 =-,44185 ~ 0820 =5,5/967 
~ 9330 =,92341 - 0830 =55 67363 
~0340 ~, 60809 - 0840 =-5, 76669 
s0320 -,69565 ,0sa0 =-5,85886 
~0360 =, 1/8588 . 9860 =5,.95012 
~0370 =-,87855 -0870 =-6, 04048 
~ 90380 ~,9/7345 - 0880 -6,.12991 
~0390 =-1,07039 - 9890 =-6, 21842 
~ 9400 =-1,16919 - 9900 =6, 30599 
~ 0410 #1, 26965 ~9910 =-6, 39264 
~0420 =1L 37162 -90920 =-6,4/834 
0430 =-1,47493 - 9930 ~§, 596309 
0440 =-1,57944 - 0940 -6, 64690 
~ 0450 =-1,68499, ~0950 -6,72977 
- 90460 =-1,79146 ~ 9960 -6,81168 
~0470 eljolJo/ 3 ~0970 =6,89263 
- 0480 =2, 00666 .9980 -~6.97264 
.0490 =2,11516 ~0990 =7,05169 
- 0500 =-2, 22412 » L000 =-7,12979 


Determination of Routh=-Hurwitz criterion for a third 
order system as sampling period, T, goes from ,001 to 0,1 
seconds, 
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